"磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出して、3軸の座標系において表現される3次元のベクトルデータとして磁気データを順次出力する3次元磁気センサを含む、複数のセンサと、前記磁気データのうち前記部品の発生する磁界の成分を表す3次元のベクトルデータであるオフセットを推定する状態変数を含む複数の状態変数を要素とするベクトルを状態ベクトルとし、前記複数のセンサの出力を要素とする観測値ベクトル及び前記状態ベクトルを用いて算出される観測残差に基づいて、前記状態ベクトルを更新するカルマンフィルタと、前記磁気データを蓄積する蓄積手段と、前記磁気データが所定数蓄積された場合に、前記所定数の前記磁気データの各々が示す3軸の座標が、第1球面の近傍に確率的に分布すると仮定して、前記第1球面の中心点を示す3軸の座標を算出する、中心点算出手段と、を備え、前記カルマンフィルタは、前記オフセットを推定する状態変数を、前記第1球面の中心点を示す3軸の座標によって更新する、ことを特徴とする、状態推定装置。磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出して、3軸の座標系において表現される3次元のベクトルデータとして磁気データを順次出力する3次元磁気センサを含む、複数のセンサと、前記磁気データのうち前記部品の発生する磁界の成分を表す3次元のベクトルデータであるオフセットを推定する状態変数を含む複数の状態変数を要素とするベクトルを状態ベクトルとし、状態遷移モデルを用いて、ある時刻における前記状態ベクトルから当該ある時刻から単位時間後の前記状態ベクトルを推定し、観測モデルを用いて、前記単位時間後の前記状態ベクトルから、推定観測値ベクトルを算出し、前記推定観測値ベクトルと前記複数のセンサの出力を要素とする観測値ベクトルとに基づいて観測残差を算出し、前記観測残差と前記単位時間後の前記状態ベクトルとに基づいて前記状態ベクトルを更新するカルマンフィルタと、前記磁気データを蓄積する蓄積手段と、前記磁気データが所定数蓄積された場合に、前記所定数の前記磁気データの各々が示す3軸の座標の分布の3次元的な広がりの程度を示す分散評価値を算出する、分布判定手段と、前記磁気データが所定数蓄積され、且つ、前記分散評価値が、分散許容値以上である場合に、前記所定数の前記磁気データの各々が示す3軸の座標が、第1球面の近傍に確率的に分布すると仮定して、前記第1球面の中心点を示す3軸の座標を算出する、中心点算出手段と、複数の前記磁気データの各々が示す3軸の座標が、第2球面と曲面とを重ね合わせることで得られる立体の表面近傍に確率的に分布すると仮定して、複数の前記磁気データに基づいて、前記立体と、前記第2球面との形状の相違の程度を示す歪評価値を算出し、前記歪評価値が歪許容値以下であるか否かを判定する、歪判定手段と、前記歪判定手段の判定結果が肯定である場合に、前記第1球面の中心点を示す3軸の座標を出力する、中心点出力手段と、を備え、前記カルマンフィルタは、前記中心点出力手段が、前記第1球面の中心点を示す3軸の座標を出力する場合に、前記オフセットを推定する状態変数を、前記第1球面の中心点を示す3軸の座標によって更新することを特徴とする、状態推定装置。複数の前記磁気データの各々で示される3軸の座標として特定されるそれぞれの位置が、前記第1球面近傍に確率的に分布すると仮定した場合の、複数の前記磁気データで示される3軸の座標により特定されるそれぞれの位置と前記第1球面との誤差を表すベクトルを、第1球面誤差ベクトルとし、3次元のベクトルを変数とし、前記第1球面誤差ベクトルの大きさを表す関数を、中心点算出関数としたとき、前記第1球面の中心点を示す3軸の座標は、前記中心点算出関数の値を最小化するときの、前記3次元のベクトルが示す3軸の座標であることを特徴とする、請求項1または2に記載の状態推定装置。磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出して、3軸の座標系において表現される3次元のベクトルデータとして磁気データを順次出力する3次元磁気センサを含む、複数のセンサと、前記磁気データを蓄積する蓄積手段と、を備える状態推定装置に用いられるオフセット更新方法であって、前記磁気データのうち前記部品の発生する磁界の成分を表す3次元のベクトルデータであるオフセットを推定する状態変数を含む複数の状態変数を要素とするベクトルを状態ベクトルとし、前記複数のセンサの出力を要素とする観測値ベクトル及び前記状態ベクトルを用いて算出される観測残差に基づいて、前記状態ベクトルを更新し、前記磁気データが所定数蓄積された場合に、前記所定数の前記磁気データの各々が示す3軸の座標が、第1球面の近傍に確率的に分布すると仮定して、前記第1球面の中心点を示す3軸の座標を算出し、前記オフセットを推定する状態変数を、前記第1球面の中心点を示す3軸の座標によって更新する、ことを特徴とする、オフセット更新方法。磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出して、3軸の座標系において表現される3次元のベクトルデータとして磁気データを順次出力する3次元磁気センサを含む、複数のセンサと、前記磁気データを蓄積する蓄積手段と、を備える状態推定装置に用いられるオフセット更新プログラムであって、前記磁気データのうち前記部品の発生する磁界の成分を表す3次元のベクトルデータであるオフセットを推定する状態変数を含む複数の状態変数を要素とするベクトルを状態ベクトルとし、前記複数のセンサの出力を要素とする観測値ベクトル及び前記状態ベクトルを用いて算出される観測残差に基づいて、前記状態ベクトルを更新する処理と、前記磁気データが所定数蓄積された場合に、前記所定数の前記磁気データの各々が示す3軸の座標が、第1球面の近傍に確率的に分布すると仮定して、前記第1球面の中心点を示す3軸の座標を算出する処理と、前記オフセットを推定する状態変数を、前記第1球面の中心点を示す3軸の座標によって更新する処理とを、コンピュータに実行させるオフセット更新プログラム。一般的に、カルマンフィルタは、動的システムの状態を表す複数の物理量の経時的な変化を推定する状態遷移モデルと、推定された動的システムの状態から、動的システムの有する複数のセンサが計測する値観測値を推定する観測モデルとを有する。そして、カルマンフィルタは、推定された観測値と、複数のセンサが実際に測定する観測値との差分観測残差を用いて、動的システムの状態を表す複数の物理量を要素とする状態ベクトルを、より真値に近い値へと更新する。しかし、観測残差は、更新前の状態ベクトルを用いて算出される値であるため、動的システムの状態が急激且つ大きく変化した場合には、状態ベクトルが真値から離れた値に更新されてしまうという問題が存在した。本発明は、上述した点に鑑みてなされたものであり、動的システムの状態が急激且つ大きく変化した場合にも、動的システムの状態を正確に推定することを解決課題とする。以下、本発明について説明する。なお、本発明の理解を容易にするために本実施形態、変形例、及び添付図面の参照符号を括弧書きにて付記するが、それにより本発明が本実施形態に限定されるものではない。上述した課題を解決するため、本発明に係る状態推定装置は、磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出して、3軸の座標系において表現される3次元のベクトルデータとして磁気データqiを順次出力する3次元磁気センサを含む、複数のセンサと、前記磁気データのうち前記部品の発生する磁界Biの成分を表す3次元のベクトルデータであるオフセットqOFFを推定する状態変数qOを含む複数の状態変数を要素とするベクトルを状態ベクトルxkとし、前記複数のセンサの出力を要素とする観測値ベクトルyk及び前記状態ベクトルを用いて算出される観測残差zkに基づいて、前記状態ベクトルを更新するカルマンフィルタと、前記磁気データを蓄積する蓄積手段と、前記磁気データが所定数N蓄積された場合に、前記所定数の前記磁気データの各々が示す3軸の座標が、第1球面Sの近傍に確率的に分布すると仮定して、前記第1球面の中心点cOを示す3軸の座標を算出する、中心点算出手段と、を備え、前記カルマンフィルタは、前記オフセットを推定する状態変数を、前記第1球面の中心点を示す3軸の座標によって更新する、ことを特徴とする。この発明によれば、状態推定装置は、中心点算出手段と、カルマンフィルタとを備え、カルマンフィルタがその演算に用いる状態ベクトルの要素のうち3次元磁気センサのオフセットを推定する状態変数を、中心点算出手段が算出する中心点の座標によって更新したうえで、カルマンフィルタの演算を行う。カルマンフィルタは、ある時刻における動的システムの状態を表す物理量状態変数を要素とする状態ベクトルを、異種の物理量を測定する複数のセンサからの出力値を要素とする観測値ベクトルと状態ベクトルとを用いて算出される観測残差に基づいて、動的システムの真の値へと近付けるように更新する演算である。具体的には、カルマンフィルタは、動的システムの状態の経時的な変化を表す状態遷移モデルを用いて、ある時刻における状態ベクトルから、当該ある時刻から単位時間経過後の状態ベクトルを推定する。次に、カルマンフィルタは、動的システムの状態の推定値から複数のセンサからの出力値を推定する観測モデルを用いて、推定された状態ベクトルから観測値ベクトルを推定する。さらに、カルマンフィルタは、推定された観測値ベクトルと、各種センサからの実際の出力実際の観測値を要素とする観測値ベクトルとの差分を演算することにより、観測残差を算出する。そして、カルマンフィルタは、観測残差を用いて、状態ベクトルを更新する。ところで、3次元磁気センサが搭載される機器は、着磁性を有する各種金属や、電気回路等、磁界を発生させる部品が備えられることが多い。この場合、3次元磁気センサが出力するベクトルデータは、地磁気を表すベクトルの他に、機器に搭載された部品が発する磁界等を表すベクトルも含む値となる。従って、地磁気の値を正確に知るためには、3次元磁気センサが出力するベクトルデータから、機器の部品が発する内部磁界を表すベクトルを取り除く補正処理が必要となる。このように、検出対象である地磁気の正確な値を得るために、補正処理において、3次元磁気センサから出力されるデータから取り除かれる内部磁界を表すベクトルを、3次元磁気センサのオフセットと呼ぶ。内部磁界は、機器の内部状態が変化しない場合には、一定の方向及び大きさを有するが、機器の内部状態が変化した場合には、方向及び大きさが変化する。機器の内部状態は、機器の利用者による機器の操作や、機器の外部の環境等に依存して、急激且つ大きな変化を生ずる場合がある。そして、機器の内部状態の変化に伴い、機器に備えられた部品が発する内部磁界すなわち、3次元磁気センサのオフセットも、急激且つ大きな変化をする。前述の通り、カルマンフィルタが行う演算には、動的システムの状態の経時的な変化を表す状態遷移モデルを用いて、ある時刻における状態変数の値から、当該ある時刻よりも単位時間経過後の状態変数の値を推定する演算が含まれる。従って、カルマンフィルタが推定する状態変数は、状態遷移モデルにおいて定式化することのできる物理量であることが好ましい。しかし、3次元磁気センサのオフセットは、機器の利用者による機器の操作等により急激且つ大きな変化をするため、この値の変化を定式化することは困難である。また、状態ベクトルを更新するために用いられる観測残差は、更新前の状態ベクトルの値に基づいて算出される。つまり、カルマンフィルタの演算により算出される更新された状態ベクトルの値は、更新前の状態ベクトルの値の影響を受けた値となる。そして、更新前の状態ベクトルの値は、それ以前の状態ベクトルの値の影響を受けた値となる。従って、カルマンフィルタの演算を繰り返す場合、カルマンフィルタが推定する状態ベクトルの値は、カルマンフィルタの演算を開始した時点から、状態ベクトルの値が推定される時点に至る全てのステップにおける状態ベクトルの値を考慮して算出される値である。よって、3次元磁気センサのオフセットが、急激且つ大きな変化をした場合には、3次元磁気センサのオフセットを推定する状態変数の更新後の値は、変化後の3次元磁気センサのオフセットから大きく乖離した値となることがある。状態ベクトルの値が、実際の物理量を正確に表した値真値から大きく乖離した値に収束している場合、その後、カルマンフィルタによる演算を繰り返しても、状態ベクトルの値は真値に近づかなくなることがあり、また、仮に状態ベクトルの値が真値に近づくことがあっても、真値に近い値に収束するまでに長い時間を要する。この場合、状態推定装置は、システムの状態を正確に推定できないため、状態推定装置が算出する状態変数に基づいて、地磁気の向きや機器の姿勢等を算出することはできない。地磁気は、磁極北に向かう水平成分と伏角方向の鉛直成分とを有する磁界であり、地面に対して一定の方向と一定の大きさとを有する一様な磁界である。よって、地面に対して機器の姿勢を変化させる場合には、機器から見た地磁気の方向も変化することになる。すなわち、機器に搭載された3次元磁気センサから見た場合、地磁気は、機器の姿勢の変化に伴い向きが変化する一定の大きさを持ったベクトルとして表される。一方、内部磁界は、機器の内部状態が変化しない場合には、機器に対して一定の方向を向き、一定の大きさを有する磁界である。よって、機器の内部状態が変化しない場合、内部磁界は、機器に搭載された3次元磁気センサから見て、機器をどのような姿勢に変化させた場合であっても、一定の方向と一定の大きさを有するベクトルとして表される。従って、3次元磁気センサの姿勢を変化させつつ取得した複数の磁気データは、内部磁界の方向及び大きさを表すベクトルの先端を中心とし、地磁気の大きさを半径とする球面の近傍に分布することになる。中心点算出手段は、複数の磁気データの示す座標の各々が、ある球面第1球面の近傍に分布すると仮定し、当該第1球面の中心点の座標を算出することにより、内部磁界、すなわち、3次元磁気センサのオフセットを算出する。中心点算出手段は、3次元磁気センサのオフセットが急激且つ大きな変化をした場合であっても、その後に取得された複数の磁気データにより正確なオフセットを算出することが可能である。本発明の状態推定装置によれば、状態ベクトルを構成する要素のうち3次元磁気センサのオフセットを推定する状態変数を、中心点算出手段により出力した中心点により上書きする。従って、3次元磁気センサのオフセットが急激且つ大きな変化をした場合であっても、変化後のオフセットを表す中心点の座標を中心点算出手段が算出し、3次元磁気センサのオフセットを推定する状態変数を中心点の座標によって上書きするため、3次元磁気センサのオフセットを推定する状態変数が真値3次元磁気センサのオフセットから乖離した状態が長期間にわたって継続することを防止することができる。すなわち、本発明に係る状態推定装置は、システムの状態を正確に推定するため、状態推定装置が算出する状態変数に基づいて、正確な地磁気の向きや、正確な機器の姿勢等を算出することが可能となる。また、上述した状態推定装置において、磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出して、3軸の座標系において表現される3次元のベクトルデータとして磁気データを順次出力する3次元磁気センサを含む、複数のセンサと、前記磁気データのうち前記部品の発生する磁界の成分を表す3次元のベクトルデータであるオフセットqOFFを推定する状態変数qOを含む複数の状態変数を要素とするベクトルを状態ベクトルxkとし、状態遷移モデルfを用いて、ある時刻における前記状態ベクトルxk-1から当該ある時刻から単位時間後の前記状態ベクトルx-kを推定し、観測モデルhを用いて、前記単位時間後の前記状態ベクトルx-kから、推定観測値ベクトルy-kを算出し、前記推定観測値ベクトルy-kと前記複数のセンサの出力を要素とする観測値ベクトルykとに基づいて観測残差zkを算出し、前記観測残差と前記単位時間後の前記状態ベクトルとに基づいて前記状態ベクトルxkを更新するカルマンフィルタと、前記磁気データを蓄積する蓄積手段と、前記磁気データが所定数N蓄積された場合に、前記所定数の前記磁気データの各々が示す3軸の座標の分布の3次元的な広がりの程度を示す分散評価値を算出する、分布判定手段と、前記磁気データが所定数蓄積され、且つ、前記分散評価値が、分散許容値λ0以上である場合に、前記所定数の前記磁気データの各々が示す3軸の座標が、第1球面Sの近傍に確率的に分布すると仮定して、前記第1球面の中心点cOを示す3軸の座標を算出する、中心点算出手段と、複数の前記磁気データの各々が示す3軸の座標が、第2球面S2と曲面SXとを重ね合わせることで得られる立体SDの表面近傍に確率的に分布すると仮定して、複数の前記磁気データに基づいて、前記立体と、前記第2球面との形状の相違の程度を示す歪評価値gDEを算出し、前記歪評価値が歪許容値δ0以下であるか否かを判定する、歪判定手段と、前記歪判定手段の判定結果が肯定である場合に、前記第1球面の中心点を示す3軸の座標を出力する、中心点出力手段と、を備え、前記カルマンフィルタは、前記中心点出力手段が、前記第1球面の中心点を示す3軸の座標を出力する場合に、前記オフセットを推定する状態変数を、前記第1球面の中心点を示す3軸の座標によって更新することを特徴とすることが好ましい。前述の通り、中心点算出手段は、複数の磁気データの示す座標が第1球面の近傍に分布すると仮定して、第1球面の中心点の座標を算出する。従って、複数の磁気データが3次元的な広がりを有さず平面的に分布する場合には、第1球面の中心点の座標を正確に特定することはできない。この発明によれば、まず、分布判定手段により複数の磁気データの3次元的な広がりの程度を示す分散評価値を算出し、その後、分散評価値が分散許容値以上である場合に、中心点算出手段が中心点の座標を算出する。従って、複数の磁気データが2次元的に分布する場合に、中心点算出手段が不正確な中心点の座標を算出することを防止することが可能となる。これにより、不正確なオフセットを用いた補正処理が施されて不正確な地磁気の方向及び大きさが算出されることを防止することが可能となる。また、不正確な中心点の座標を算出する演算処理を未然に防止することが可能となるため、状態推定装置の低消費電力化が可能となった。また、この発明によれば、歪判定手段により、複数の磁気データの座標を表面近傍に有する立体と、第2球面との形状の相違の程度を、歪評価値に基づいて判定する。歪評価値が閾値よりも大きな値を有する場合には、立体の形状が第2球面とは異なる歪んだ形状を有するため、立体の表面近傍に分布する複数の磁気データの示す座標は、第1球面の近傍に分布すると看做すことはできない。中心点算出手段は、第1球面の近傍に複数の磁気データが存在していることを前提として、3次元磁気センサのオフセットの候補である中心点の座標を算出するため、複数の磁気データを近傍に有さない第1球面の中心点の座標を示すベクトルは、3次元磁気センサのオフセットとしての意味を有さない。本発明に係る状態推定装置は、歪判定手段を備えるため、3次元磁気センサのオフセットと看做すことのできない不適切な中心点の座標を示すベクトルにより、3次元磁気センサのオフセットを推定する状態変数を上書きすることを防止することが可能となった。このように、本発明に係る状態推定装置は、分布判定手段及び歪判定手段を備えるため、3次元磁気センサのオフセットと看做すことのできる真値に近い中心点の座標が得られた場合にのみ、当該真値に近い中心点の座標を示すベクトルにより、3次元磁気センサのオフセットを推定する状態変数を上書きする。これにより、本発明に係る状態推定装置は、システムの状態を正確に推定することができ、状態推定装置が算出する状態変数に基づいて、正確な地磁気の向きや、正確な機器の姿勢等を算出することが可能となる。また、本発明の具体的な態様として、上述した状態推定装置は、前記歪判定手段の判定結果が否定である場合に、利用者に前記機器の位置を変化させずに前記機器の姿勢を変化させるように促す手段を備えてもよい。この発明によれば、上述した状態推定装置は、歪判定手段において歪評価値が歪許容値よりも大きいと判定される場合、すなわち、3次元磁気センサから出力される複数の磁気データの示す座標が立体の表面近傍に分布するときに立体の形状が球面から大きな歪を有する場合に、利用者に対して機器の位置を変化させずに機器の姿勢を変化させるように促す手段を備える。3次元磁気センサは、測定対象である地磁気の他に、機器に備えられた部品の発する内部磁界と、外部の物体の発する不均一な外部磁界とを検出する場合がある。そして、この不均一な外部磁界の影響が大きい場合には、複数の磁気データの示す座標は、球面とは大きく異なる歪んだ形状を有する立体の表面近傍に分布することになり、中心点算出手段により算出された中心点の座標を示すベクトルを3次元磁気センサのオフセットとして採用することができなくなる。しかし、このような不均一な外部磁界も、外部磁界を発生させる外部の物体と3次元磁気センサとの相互位置関係が変化しない場合には、一定の大きさの磁界にすぎない。つまり、3次元磁気センサの位置を変化させずに姿勢を変化させる場合、外部磁界は、3次元磁気センサによって、一定の大きさを有し方向のみを変化させる磁界として検出される。この場合、3次元磁気センサから出力される複数の磁気データで示される座標は、ある球面近傍に分布することになる。そして、この球面の中心点の座標は、機器に備えられた部品の発する磁界内部磁界の成分を表すベクトルが示す座標とほぼ同一の座標となる。前述の通り、第1球面は、複数の磁気データの示す座標が第1球面の近傍に分布することを仮定して定められる。よって、複数の磁気データの示す座標が球面の形状となるように分布する場合、複数の磁気データの示す座標は、中心点算出手段によって算出された中心点を中心とする第1球面近傍にも分布すると看做すことができる。従って、第1球面の中心点は内部磁界の成分を表すベクトルが示す座標とほぼ同一の座標となるため、第1球面の中心点の座標を示すベクトルを3次元磁気センサのオフセットとして採用することが可能となる。このように、上述した状態推定装置は、利用者に機器の位置を変化させずに機器の姿勢を変化させるように促す手段を備えることにより、不均一な外部磁界が存在する場所であっても、適切なオフセットを算出することができるような動作を利用者に促すことが可能となる。また、本発明の具体的な態様として、上述した状態推定装置において、前記分散評価値は、複数の前記磁気データの各々が示す3軸の座標の分散を表す共分散行列Aの最小固有値λ3であることとしてもよい。この発明によれば、複数の磁気データの分散を表す共分散行列の最小固有値の値を、分散評価値として用いている。磁気データは3軸の座標を示すので、複数の磁気データの分散を表す共分散行列は、3行3列の行列となる。この共分散行列から3つの固有ベクトルと3つの固有値が得られる。3つの固有ベクトルの一つは、複数の磁気データが最も大きく分布する方向を表している。この固有ベクトルに対応する固有値は最大固有値となり、逆に、複数の磁気データが最も小さく分布する方向の固有ベクトルに対応する固有値が最小固有値となる。複数の磁気データが理想的に2次元的に分布する場合、最小固有値の値は0に限りなく近い値を有し、複数の磁気データが3次元的に分布する場合、最小固有値の値は、複数の磁気データの分布の3次元的な広がりの程度に従って大きな値を有する。従って、共分散行列の最小固有値を、分散評価値として用いることにより、複数の磁気データの分布の3次元的な広がりの程度を的確に把握することが可能となる。また、本発明の具体的な態様として、上述した状態推定装置において、前記磁気データを前記第1球面の中心点を原点とする座標系で表した3軸のベクトルqi-cOと、当該3軸のベクトルを対称行列である歪評価行列Eにより変換したベクトルEとの内積を、複数の前記磁気データの各々について算出し、複数の算出結果を要素とするベクトルを歪誤差ベクトルkとし、複数の前記磁気データで示される3軸の座標により特定されるそれぞれの位置と前記第2球面S2との誤差を表すベクトルを、第2球面誤差ベクトルδS2とし、前記歪誤差ベクトルと、前記第2球面誤差ベクトルとの和を立体誤差ベクトルδSDとし、前記歪評価行列Eの各々の成分e11~e33と、前記第2球面の中心点cO2を示す3軸の座標とを変数とし、前記立体誤差ベクトルの大きさを表す関数を歪評価関数fSDとしたとき、前記歪評価値は、前記歪評価関数を最小化するときの、前記歪評価行列のノルムであることとしてもよい。この発明によれば、立体誤差ベクトルの大きさを最小化することで、歪評価行列と、第2球面の中心点の座標とを算出する。そして、歪評価値は、歪評価行列のノルムとして算出される。歪評価行列は、3軸のベクトルを座標変換するための対称行列であるため、3行3列の行列である。従って、歪評価値は3つの固有値を有し、3つの固有値のうち、絶対値の大きさが最大となる固有値の絶対値が、歪評価値となる。立体誤差ベクトルは、第2球面誤差ベクトルと、歪誤差ベクトルとの和で与えられる。第2球面誤差ベクトルは、単に、複数の磁気データで示される座標と第2球面との誤差を表すためのベクトルである。従って、複数の磁気データで示される座標と第2球面との誤差を最小化するように第2球面を定めた場合、球面誤差ベクトルにより表される複数の誤差は、全体として対称性を有し、方向依存性の無い、ホワイトノイズとなる。一方、歪誤差ベクトルは、第1球面の中心点から見た磁気データの示す座標を表すベクトルと、当該ベクトルを歪評価行列により変換したベクトルとの内積を、複数の磁気データの各々について並べたベクトルである。換言すれば、歪誤差ベクトルの各要素は、対称行列である歪評価行列を係数行列とし、中心点から見た磁気データの座標を表すベクトルの3つの要素の各々を変数とする、3変数二次形式となる。すなわち、歪誤差ベクトルは、複数の磁気データで示される座標と第2球面との誤差の各々を、二次形式で表される同一の関数に基づく曲面上に存在するという制約の下で表現するベクトルである。このような歪誤差ベクトルを用いて、複数の磁気データと第2球面との誤差を表現することにより、ホワイトノイズとは別に、二次形式の関数に基づく曲面で表される誤差すなわち、球面からの歪を表現することが可能となる。球面誤差ベクトルと歪誤差ベクトルとを加算して得られる立体誤差ベクトルは、複数の磁気データを表面近傍に有し、球面と曲面とを重ね合わせた球面とは異なる立体を表現する。そして、立体誤差ベクトルのうち、歪誤差ベクトルの大きさ評価することで、立体の形状が、第2球面の形状からどの程度の差異歪を有するかつまり、ホワイトノイズとは異なる球面からの歪による誤差の大きさを、評価することが可能となる。その結果、立体と第2球面とが実質的に同一の形状を有するほど近似していると評価される場合には、複数の磁気データの示す座標が第2球面近傍にも分布すると看做すことができる。前述の通り、複数の磁気データの示す座標が、球面の形状となるように分布する場合、複数の磁気データの示す座標は、第1球面近傍にも分布すると看做すことができる。従って、立体と第2球面とが実質的に同一の形状を有するほど近似していると評価される場合には、第1球面の中心点の座標を示すベクトルを3次元磁気センサのオフセットとして採用することができる。一方、立体が第2球面とは大きく異なる歪んだ形状を有すると評価される場合には、複数の磁気データで示される座標が第2球面近傍に分布すると看做すことはできないため、第1球面の中心点の座標を示すベクトルを3次元磁気センサのオフセットとして採用することを防止する必要がある。すなわち、上述した状態推定装置は、中心点算出手段によって算出された中心点の座標のうち、不均一な外部磁界の影響を受けていない複数の磁気データより算出された、適切な中心点の座標を示すベクトルのみを、3次元磁気センサのオフセットとして採用することが可能である。この結果、上述した状態推定装置は、システムの状態を正確に推定するため、状態推定装置が算出する状態変数に基づいて、正確な地磁気の向きや、正確な機器の姿勢等を算出することが可能となる。また、上述した状態推定装置において、複数の前記磁気データの各々で示される3軸の座標として特定されるそれぞれの位置が、前記第1球面近傍に確率的に分布すると仮定した場合の、複数の前記磁気データで示される3軸の座標により特定されるそれぞれの位置と前記第1球面との誤差を表すベクトルを、第1球面誤差ベクトルδSとし、3次元のベクトルcを変数とし、前記第1球面誤差ベクトルδSの大きさを表す関数を、中心点算出関数fScとしたとき、前記第1球面の中心点を示す3軸の座標は、前記中心点算出関数の値を最小化するときの、前記3次元のベクトルが示す3軸の座標であることを特徴とすることが好ましい。この発明によれば、簡易な計算により、3次元磁気センサのオフセットの候補となる座標を算出することが可能となる。また、本発明の具体的な態様として、上述した状態推定装置において、前記歪評価関数は、下記のfSDで表すこともできる。この発明によれば、歪評価関数を最小化するときの歪評価行列のノルムである歪評価値を簡単な計算により算出することができる。従って、上述した状態推定装置は、処理速度の向上、及び低消費電力化が可能になるという利点を有する。また、本発明の具体的な態様として、上述した状態推定装置において、前記中心点算出関数は、下記のfScで表すこともできる。この発明によれば、簡単な計算により中心点を算出することができる。従って、上述した状態推定装置は、処理速度の向上、及び低消費電力化が可能になるという利点を有する。次に、本発明に係る状態推定装置に用いられるオフセット更新方法は、磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出して、3軸の座標系において表現される3次元のベクトルデータとして磁気データを順次出力する3次元磁気センサを含む、複数のセンサと、前記磁気データを蓄積する蓄積手段と、を備える状態推定装置に用いられるオフセット更新方法であって、前記磁気データのうち前記部品の発生する磁界の成分を表す3次元のベクトルデータであるオフセットを推定する状態変数を含む複数の状態変数を要素とするベクトルを状態ベクトルとし、前記複数のセンサの出力を要素とする観測値ベクトル及び前記状態ベクトルを用いて算出される観測残差に基づいて、前記状態ベクトルを更新し、前記磁気データが所定数蓄積された場合に、前記所定数の前記磁気データの各々が示す3軸の座標が、第1球面の近傍に確率的に分布すると仮定して、前記第1球面の中心点を示す3軸の座標を算出し、前記オフセットを推定する状態変数を、前記第1球面の中心点を示す3軸の座標によって更新する、ことを特徴とする。この発明によれば、複数の磁気データに基づいて算出された中心点の座標により、状態ベクトルを構成する要素のうち3次元磁気センサのオフセットを推定する状態変数を上書きする。従って、3次元磁気センサのオフセットが急激且つ大きな変化をした場合であっても、3次元磁気センサのオフセットを推定する状態変数が、真値3次元磁気センサのオフセットから乖離することを防止することが可能となる。すなわち、状態推定装置は、システムの状態を正確に推定するため、状態推定装置が算出する状態変数に基づいて、正確な地磁気の向きや、正確な機器の姿勢等を算出することが可能となる。次に、本発明に係る状態推定装置に用いられるオフセット更新プログラムは、磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出して、3軸の座標系において表現される3次元のベクトルデータとして磁気データを順次出力する3次元磁気センサを含む、複数のセンサと、前記磁気データを蓄積する蓄積手段と、を備える状態推定装置に用いられるオフセット更新プログラムであって、前記磁気データのうち前記部品の発生する磁界の成分を表す3次元のベクトルデータであるオフセットを推定する状態変数を含む複数の状態変数を要素とするベクトルを状態ベクトルとし、前記複数のセンサの出力を要素とする観測値ベクトル及び前記状態ベクトルを用いて算出される観測残差に基づいて、前記状態ベクトルを更新する処理と、前記磁気データが所定数蓄積された場合に、前記所定数の前記磁気データの各々が示す3軸の座標が、第1球面の近傍に確率的に分布すると仮定して、前記第1球面の中心点を示す3軸の座標を算出する処理と、前記オフセットを推定する状態変数を、前記第1球面の中心点を示す3軸の座標によって更新する処理とを、コンピュータに実行させることを特徴とする。この発明によれば、複数の磁気データに基づいて算出された中心点の座標により、状態ベクトルを構成する要素のうち3次元磁気センサのオフセットを推定する状態変数を上書きする。従って、3次元磁気センサのオフセットが急激且つ大きな変化をした場合であっても、3次元磁気センサのオフセットを推定する状態変数が、真値3次元磁気センサのオフセットから乖離することを防止することが可能となる。すなわち、状態推定装置は、システムの状態を正確に推定するため、状態推定装置が算出する状態変数に基づいて、正確な地磁気の向きや、正確な機器の姿勢等を算出することが可能となる。<A.実施形態>本発明の実施の形態について図面を参照して説明する。[1.機器構成及びソフトウェア構成]図1は、本発明の実施形態に係る携帯機器1のブロック図であり、図2は携帯機器1の外観を示す斜視図である。携帯機器1は、画面の姿勢に応じて地図などの画像を回転させることによって、画像によって表されている方位を現実空間の方位に追従させる機能を備える。この機能は、各種センサの出力に基づいてカルマンフィルタの演算を行い、携帯機器1の姿勢と、携帯機器1から見た地磁気の方向とを推定することによって実現される。携帯機器1は、各種の構成要素とバスを介して接続され装置全体を制御するCPU10、CPU10の作業領域として機能するRAM20、各種のプログラムやデータを記憶したROM30、通信を実行する通信部40、画像を表示する表示部50、及びGPS部60を備える。GPS部60は、GPS衛星からの信号を受信して携帯機器1の位置情報緯度、経度を生成するものである。また、携帯機器1は、地磁気を検出して磁気データを出力する3次元磁気センサ70、加速度を検出して加速度データを出力する3次元加速度センサ80、及び角速度を検出して角速度データを出力する3次元角速度センサ90を備える。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は、磁極北に向かう水平成分と伏角方向の鉛直成分とを有する地磁気を表すベクトルの他に、携帯機器1に搭載された部品が発する磁界等を表すベクトルも含む値となる。従って、地磁気の値を正確に知るためには、3次元磁気センサが出力するベクトルデータ磁気データqから、携帯機器1の部品が発する内部磁界を表すベクトルを取り除く補正処理が必要となる。このように、検出対象である地磁気の正確な値を得るために、補正処理において、3次元磁気センサ70から出力されるデータから取り除かれる内部磁界を表すベクトルを磁気センサのオフセットqOFFと呼ぶ。また、3次元磁気センサ70は、図3に示すように、地磁気Bg及び内部磁界Biの他に、携帯機器1の外部に存在する物体2が出す外部磁界Bxを検出する場合がある。詳細は後述するが、外部磁界Bxは、不均一な磁界である。従って、外部磁界Bxの成分を算出して、これを3次元磁気センサ70の出力から取り除いて、正確な地磁気Bgの値を算出することは困難である。一方、外部磁界Bxが存在しない場合、または外部磁界Bxの影響の小さな場合には、3次元磁気センサ70の出力から、内部磁界Biをオフセットとして取り除く処理を施すことにより、正確な地磁気Bgを算出することが可能である。ここで、説明の便宜上、図3に示すような地上座標系ΣG及びセンサ座標系ΣSを導入する。地上座標系ΣGは、地上に固定された座標系であり、地上の任意の一点を原点とし、互いに直交する3つの方向、例えば、水平東、水平北、及び鉛直上向きを、それぞれx軸、y軸、及びz軸とする座標系である。センサ座標系ΣSは、携帯機器1に固定された座標系である。3次元磁気センサ70が出力する磁気データは、センサ座標系ΣSのx軸、y軸、z軸上にプロットされ、センサ座標系ΣSのベクトルデータとして表現される。また、地上座標系ΣGから見たセンサ座標系ΣSの原点の位置及び姿勢つまり、地上座標系ΣGから見た携帯機器1の位置及び姿勢を、それぞれ、位置PS及び姿勢μと表現する。なお、図3に記載された各ベクトルの左上に付された添字Gは、当該ベクトルが地上座標系ΣGにおいて表現されたベクトルであることを意味する。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に格納されている状態推定プログラムを実行することによって、携帯機器1の姿勢や磁気センサのオフセット等の物理量状態変数を推定する。すなわち、携帯機器1は、CPU10が状態推定プログラムを実行することにより、状態推定装置として機能する。図4は、状態推定装置のうち、CPU10が状態推定プログラムを実行することにより実現される機能を表す機能ブロック図である。図4に示すように、状態推定装置は、状態推定部100と、中心点導出部200とを備える。状態推定部100は、カルマンフィルタ演算部120と、初期値生成部140とを備える。カルマンフィルタ演算部120は、磁気データq、加速度データa、及び角速度データgを要素とする観測値ベクトルyを用いてカルマンフィルタの演算を実行することで、状態ベクトルxを周期的に更新し、更新後の状態ベクトルxを出力する。初期値生成部140は、各種設定情報に基づいて初期値INIを生成し、生成した初期値INIを出力する。中心点導出部200は、3次元磁気センサ70より供給される磁気データqを複数個蓄積したうえで、複数個の磁気データq1~qNに基づいて中心点第1球面の中心点cOの座標を算出し、算出した中心点cOの座標を出力するNは精度よく中心点cOの座標を導出するために必要な磁気データの規定測定回数を表す5以上の自然数。詳細は後述するが、中心点cOは、センサ座標系ΣS上の座標を示す点であり、中心点cOの座標を示すベクトルは、磁気センサのオフセットqOFFの候補となる。[2.カルマンフィルタについて]図5は、状態推定部100の動作を表すフローチャートである。当該フローチャートに示される状態推定処理は、CPU10が状態推定プログラムを実行することにより開始される。ステップS101において、CPU10は、初期値生成処理を実行する。具体的には、CPU10は、ROM30に記録された各種設定情報を参照し、初期値INIを生成する。すなわち、CPU10は、ステップS101の初期値生成処理を実行することにより、初期値生成部140として機能する。ステップS102において、CPU10は、中心点導出有無判定処理を実行する。具体的には、CPU10は、後述する中心点導出処理図7参照において中心点cOの座標が生成されたか否かを判定するより具体的には、CPU10は、中心点cOの座標が、後述するステップS206の中心点出力処理から引き渡されたか否かを判定する。そして、中心点導出処理において中心点cOの座標が生成された場合、CPU10は、処理をステップS103に進める。一方、中心点導出処理において中心点cOの座標が生成されていない場合、CPU10は、処理をステップS104に進める。ステップS103において、CPU10は、状態変数上書処理を実行する。詳細は後述するが、後述する状態ベクトル更新処理ステップS104において行われるカルマンフィルタの演算の対象である状態ベクトルxは、複数の状態変数を要素とするベクトルであり、複数の要素には、3次元のベクトルである磁気センサのオフセット推定値qOが含まれる。状態変数上書処理において、CPU10は、中心点導出処理によって得られた中心点cOを表す3軸の座標の各成分によって、3次元のベクトルである磁気センサのオフセット推定値qOを構成する各要素を上書きする。ステップS104において、CPU10は、状態ベクトル更新処理を実行する。具体的には、CPU10は、観測値ベクトルyを用いてカルマンフィルタの演算を実行することで、状態ベクトルxを周期的に更新し、更新後の状態ベクトルxを算出する。ステップS105において、CPU10は、状態ベクトル出力処理を実行する。具体的には、CPU10は、状態ベクトル更新処理で算出された更新後の状態ベクトルxを出力する。CPU10は、以上に示したステップS102~S105の処理を実行することにより、カルマンフィルタ演算部120として機能する。そして、携帯機器1は、CPU10がステップS102~S105の処理を実行することにより生成される更新後の状態ベクトルxkつまり、複数の状態変数を用いて、地磁気Bgの向きや、携帯機器1の姿勢等を算出することができる。次に、CPU10は、状態推定プログラムを終了する終了条件を充足するか否かを判定するステップS106。終了条件は、携帯機器1の仕様等に基づいて適宜定めればよい。例えば、携帯機器1の電源がオフ状態になることを終了条件としてもよい。終了条件を充足する場合、CPU10は、当該フローチャートに示す状態推定処理を終了する。一方、終了条件を充足しない場合、CPU10は、処理をステップS102に進める。一般的に、カルマンフィルタは、動的システムの状態を表す複数の物理量の経時的な変化を推定する状態遷移モデルと、動的システムの有する複数のセンサが計測する値観測値を推定する観測モデルと、を有する。そして、カルマンフィルタは、状態遷移モデルを用いて、ある時刻における動的システムの状態を表す複数の物理量状態変数を要素とする状態ベクトルから、単位時間が経過した後の状態ベクトルを推定する。次に、カルマンフィルタは、推定された状態ベクトルに基づいて、動的システムの状態を測定する複数のセンサの出力値を要素とする観測値ベクトルの値を推定する。さらに、推定された観測値ベクトルと、実際のセンサの出力値を要素とする観測値ベクトルとの差分として算出される観測残差とに基づいて、推定された状態ベクトルを、実際の値真値に近い値へと更新し、更新後の状態ベクトルを出力する。カルマンフィルタは、以上のような演算を繰り返すことにより、状態ベクトルを構成する複数の状態変数の各々を、実際の物理量に出来るだけ近く実際の物理量を正確に表した値真値へと近付ける。本実施形態におけるカルマンフィルタは、状態変数として、携帯機器1の姿勢μ、地磁気の強さr、地磁気の伏角φ、携帯機器1の角速度ω、角速度センサ91~93のオフセット推定値gO、及び磁気センサ71~73のオフセット推定値qOを採用し、これらを推定の対象とする。なお、状態変数に、角速度センサのオフセット推定値gOを含ませているのは、角速度データgに、角速度センサセンサ91~93のオフセットが重畳するからである。これらの状態変数を要素とする時刻T=kにおける状態ベクトルxkは、以下の式1で表される。なお、各状態変数の右下に付された添え字kは、当該状態変数が時刻T=kにおける値であることを表す。本実施形態では、姿勢μを、クォータニオンを用いて表現する。クォータニオンとは、物体の姿勢回転状態を表す4次元数である。具体的には、携帯機器1に固定されたセンサ座標系ΣSの各軸と、地上座標系ΣGの各軸とが一致する姿勢μを基準姿勢と定義し、基準姿勢をμ=0、0、0、1と表現する。このとき、携帯機器1の任意の姿勢μは、基準姿勢に対して単位ベクトルαの方向を回転軸として角度ψだけ回転した姿勢として表現できる。回転後の姿勢μは、式2で与えられる。本実施形態におけるカルマンフィルタの観測対象は、3軸の磁気センサ71~73から出力される磁気データq、3軸の加速度センサ81~83から出力される加速度データa、及び3軸の角速度センサ91~93から出力される角速度データgである。このとき、時刻T=kにおける観測値ベクトルykは式3で与えられる。前述のとおり、初期値生成部140は、時刻T=0における状態変数を要素とする初期値INIを算出する初期値生成処理。初期値INIは、姿勢μの初期値μ0、地磁気の強さrの初期値r0、地磁気の伏角φの初期値φ0、角速度ωの初期値ω0、角速度センサのオフセット推定値gOの初期値gO0、及び磁気センサのオフセット推定値qOの初期値qO0を要素として含む。初期値INIは、カルマンフィルタの演算によって状態変数がなるべく早く正確な値に収束するような値に適宜設定すればよいが、例えば、以下のような値に設定してもよい。地磁気の強さrの初期値r0、及び地磁気の伏角φの初期値φ0は、例えば、GPS部60から供給される位置情報に基づいて生成することができる。これは、地球上の位置が特定できれば、その位置における地磁気の強さ及び伏角を知ることができるからである。具体的には、ROM30には、位置情報と地磁気の強さ及び伏角とを対応づけて記憶したルックアップテーブルLUT1が格納される。初期値生成部140は、ルックアップテーブルLUT1を参照して、位置情報に対応する地磁気の強さ及び伏角を、地磁気の強さrの初期値r0及び地磁気の伏角φの初期値φ0として設定する。なお、携帯機器1が衛星からの電波の届かない場所例えば、地下街にある場合には、GPS部60から位置情報が得られない。そのような場合には、携帯機器1が使われ得る地域の典型的な値を採用すればよい。その値もルックアップテーブルLUT1に格納されている。初期値生成部140は、位置情報の取得が不能な場合には、典型的な値をルックアップテーブルLUT1から読み出す。例えば、日本で販売された携帯機器1については、明石市の地磁気の強さ及び伏角に基づいて、地磁気の強さrの初期値r0、及び地磁気の伏角φの初期値φ0を算出する。角速度ωの初期値ω0は、例えば、携帯機器1が静止していることを仮定して、0に設定する。角速度センサのオフセット推定値gOの初期値gO0は、通常0近辺に調整されているため、0に設定する。姿勢μの初期値μ0に関しては、例えば、携帯機器1が一定方向に向いて静止していることを仮定して、実際の初期姿勢とのずれを小さくするような値を設定する。磁気センサのオフセット推定値qOの初期値qO0は、例えば、時刻T=0の磁気センサの観測値q0、姿勢μの初期値μ0、携帯機器1を使用する地域の地磁気ベクトルqtypical、及び、式5に示す行列Bμを用いて、以下に示す式4で与えられる値に設定する。ここで、行列Bμは、携帯機器1が姿勢μである場合に、地上座標系ΣGにおいて表現されたベクトルを、センサ座標系ΣSで表現するための座標変換行列である。なお、ROM30には、位置情報と地磁気ベクトルqtypicalとを対応づけて記憶したルックアップテーブルLUT2が記憶されている。初期値生成部140は、GPS部60で生成される位置情報に基づいてルックアップテーブルLUT2を参照して地磁気ベクトルqtypicalを取得し、式4を演算することによって、磁気センサのオフセット推定値qOの初期値qO0を得る。初期値生成部140は、以上のようにして生成された姿勢μの初期値μ0、地磁気の強さrの初期値r0、地磁気の伏角φの初期値φ0、角速度ωの初期値ω0、角速度センサのオフセット推定値gOの初期値gO0、及び磁気センサのオフセット推定値qOの初期値qO0を要素とするベクトルである初期値INIを生成し、これを出力する。[3.カルマンフィルタによる演算]次に、ステップS104において行われる、状態ベクトル更新処理の内容、すなわち、本実施形態に係るカルマンフィルタによる演算の内容について説明する。カルマンフィルタは、動的システムの状態の経時的な変化を表す状態遷移モデルを用いて、ある時刻時刻T=k-1のシステムの状態を示す状態ベクトルxk-1から単位時間が経過した後時刻T=kの状態ベクトルxkを推定し、この推定値を、推定された状態ベクトルx-kとして出力する。そして、カルマンフィルタは、各種センサからの出力を要素とする観測値ベクトルykと状態ベクトルxkとの関係を表す観測モデルを用いて、状態遷移モデルにより推定された状態ベクトルx-kに基づいて観測値ベクトルykを推定し、この推定値を、推定された観測値ベクトル推定観測値ベクトルy-kとして出力する。なお、詳細は後述するが、本実施形態ではカルマンフィルタとして、非線形カルマンフィルタであるシグマポイントカルマンフィルタを用いる。シグマポイントカルマンフィルタとは、状態ベクトルxk-1の周囲に複数のシグマポイントχk-1を設定し、これらの複数のシグマポイントχk-1を状態遷移モデルに適用することで、単位時間経過後のシグマポイントχ-kを推定し、推定された複数のシグマポイントχ-kの平均を算出することにより、推定された状態ベクトルx-kを求める演算である。従って、推定された観測値ベクトルy-kは、厳密には、推定された状態ベクトルx-kの周辺に存在する複数のシグマポイントχ-kに基づいて算出される。次に、カルマンフィルタは、推定された観測値ベクトルy-kと、実際の観測値を要素とする観測値ベクトルykとの差分を観測残差zkとして算出し、観測残差zkに基づいてカルマンゲインKkを算出する。さらに、カルマンフィルタは、観測残差zkとカルマンゲインKkとを用いて、推定された状態ベクトルx-kを更新し、更新後の状態ベクトルxkを算出する。状態ベクトル更新処理は、以上のような状態ベクトルxkの更新を繰り返すカルマンフィルタの演算により、状態ベクトルxkを実際の物理量を正確に表した値真値に近い正確な値へと近付ける処理である。[3.1.カルマンフィルタによる演算の概要]本実施形態における、カルマンフィルタの状態遷移モデルは以下に示す式6で与えられ、観測モデルは以下に示す式7で与えられる。なお、本実施形態では、式6に現れる関数f及び式7に現れる関数hは、非線形の関数である。ここで、状態ベクトルxkはa次元のベクトルであり、観測値ベクトルykはb次元のベクトルである。なお、本実施形態では、a=15であり、b=9である。また、式6に現れるプロセスノイズwk及び式7に現れる観測ノイズvkは0を中心とするガウスノイズである。式6は、時刻T=kにおける状態ベクトルxkが、時刻T=k-1における状態ベクトルxk-1を関数fに代入して得られた値と、時刻T=k-1におけるプロセスノイズwk-1とを加算することにより推定されることを示している。また、式7は、時刻T=kにおける観測値ベクトルykが、時刻T=kにおける状態ベクトルxkを関数hに代入して得られる値と、時刻T=kにおける観測ノイズvkとを加算することにより推定されることを示している。なお、プロセスノイズwkの共分散をQk、観測ノイズvkの共分散をRkと表す。時刻T=kにおける観測残差zkは、実際の観測値ベクトルykと、推定された観測値ベクトルy-kとに基づいて定められるベクトルであり、以下に示す式8で表される。カルマンフィルタは、観測残差zkと、式9に示すカルマンゲインKkとを用いて、以下の式10に示すように、推定された状態ベクトルx-kを更新し、更新後の状態ベクトルxkを算出する。また、カルマンフィルタは、以下の式11に示すように、状態ベクトルxkの推定誤差の共分散Pkを更新する。観測モデルを用いて推定された観測値ベクトルy-kは、式7に示すように、状態遷移モデルを用いて推定された状態ベクトルx-kを、観測モデルに適用することで算出される理論値である。よって、観測モデルを用いて推定された観測値ベクトルy-kと実際のセンサ出力に基づく観測値ベクトルykとの差分である観測残差zkは、推定された状態ベクトルx-kと実際の物理量を正確に表した値真値との近似度を示す値である。式10では、このような観測残差zkを用いて、推定された状態ベクトルx-kを、真値に近い更新後の状態ベクトルxkへと更新する。[3.2.シグマポイントカルマンフィルタによる演算]なお、本実施形態では、カルマンフィルタカルマンフィルタ演算部120は、アンセンテッド変換を用いたシグマポイントカルマンフィルタにより構成される。以下では、図6に示す、カルマンフィルタ演算部120の機能ブロック図を参照しつつ、シグマポイントカルマンフィルタによる状態ベクトルの更新について、具体的に説明する。遅延部121は、加算器129から出力される更新後の状態ベクトルxkを、単位時間時刻T=k-1から時刻T=kに相当する時間だけ遅延させることで、状態ベクトルxk-1を生成し、これを、シグマポイント生成部122に対して出力する。なお、初回の演算時刻T=0では、遅延部121は、初期値生成部140が出力する初期値INIを用いて、状態ベクトルxk-1を生成する。また、中心点導出部200が中心点cOを出力する場合、遅延部121は、状態ベクトルxk-1のうち、磁気センサのオフセット推定値qOに対応する要素を、中心点導出部200より出力される中心点cOにより上書きしたうえで、状態ベクトルxk-1を生成するステップS103の状態変数上書処理。シグマポイント生成部122は、a行a列の共分散行列Pk-1及びQk-1を用いて、2a1個のシグマポイントを生成する。具体的には、まず、複数のシグマポイントの平均に対して、平均のまわりのシグマボイントの広がりを表すスケーリングパラメータλを用いて、式12及び式13に示すベクトルσkを定義する。このとき、シグマポイント生成部122は、ベクトルσk-1と状態ベクトルxk-1とに基づいて、式14及び式15で示される2a1個のシグマポイントχk-1を生成する。状態遷移モデル部123は、式16に示すように、時刻T=k-1における2a1個のシグマポイントχk-1jの各々を、状態遷移モデルに適用することにより、時刻T=kにおける2a1個のシグマポイントの推定値χ-kjを算出する。次に、平均算出部124は、式17に示すように、時刻T=kにおける2a1個のシグマポイントの推定値χ-kの平均値を計算することで、推定された状態ベクトルx-kを算出する。また、平均算出部124は、式18に示す、推定された状態ベクトルx-kの誤差の共分散P-kを算出する。一方、観測モデル部125は、式19に示すように、時刻T=kにおける2a1個のシグマポイントχkjの各々を、観測モデルに適用することにより、2a1個の推定観測値γkjを算出する。そして、平均処理部126は、式20に示すように、2a1個の推定観測値γkjの平均を演算することにより、推定された観測値ベクトルy-kを算出する。次に、平均処理部126は、式21に示す観測残差の共分散Pyykを算出する。減算器127は、式8に示したように、観測値ベクトルykと、推定された観測値ベクトルy-kとの差分として、観測残差zkを算出する。カルマンゲイン付与部128は、式22に示す相互共分散行列Pxykを算出する。そして、カルマンゲイン付与部128は、式9に示したように、残差の共分散Pyykと、相互共分散行列Pxykとに基づいて、カルマンゲインKkを算出し、式10の右辺第2項Kkzkの演算を実行する。また、カルマンゲイン付与部128は、式11に示したように、状態ベクトルxkの推定誤差の共分散Pkを、P-kからPkに更新する。加算器129は、式10に示したように、推定された状態ベクトルx-kと、カルマンゲイン付与部128から出力される式10の右辺第2項Kkzkとを加算することにより、更新後の状態ベクトルxkを算出する。[3.3.状態遷移モデル]次に、状態遷移モデル部123の演算で用いる状態遷移モデルについて説明する。本実施形態に係る状態遷移モデルにおいて、状態ベクトルxkを構成する状態変数のうち、姿勢μについての状態遷移は、式23により定義される。式23は、時刻T=k-1における姿勢μk-1から、単位時間経過後の時刻T=kにおける姿勢μkを推定する演算を表す。ここで、μーkは、時刻kにおける推定された状態ベクトルx-kのうち、姿勢μを表す状態変数に対応する要素である。なお、式23の右辺の演算子Ωは、式24により定義される。ここで、I3×3は3行3列の単位行列を表す。3次元ベクトルl=l1,l2,l3に対して、演算子[l×]は、式25で定義される。また、測定時間間隔時刻T=k-1から時刻T=kまでの間隔をΔtで表し、時刻T=kにおける更新後の状態ベクトルxkのうち角速度を表す状態変数に対応する要素をωkで表したとき、演算子Ωを構成する成分Ψkは、式26で定義される。姿勢μは、クォータニオンで表現されるため、正規化条件μ=1が満たされていないといけないが、シグマポイントの平均を求めると正規化条件が満たされなくなる。そこで、姿勢μに対して何らかの演算が行われるときには、演算後の結果をそのベクトル自身の大きさで正規化する。なお、より厳密に正規化条件を保つため、状態ベクトルxkのうち姿勢μkについてはMRPsを用いて前時刻との差分情報だけに限定し、カルマンフィルタの外部にある姿勢情報をカルマンフィルタから得られる差分情報に基づいて更新してもよい。地磁気の強さr及び地磁気の伏角φは、変化を予測することが難しい。そこで、本実施形態に係る状態遷移モデルでは、便宜上、時刻T=kにおける地磁気の強さrk及び伏角φkと、時刻T=k-1における地磁気の強さrk-1及び伏角φk-1とは等しい値であると仮定する。同様に、角速度センサのオフセット推定値gOは、変化を予測することが難しい。そこで、本実施形態に係る状態遷移モデルでは、便宜上、時刻T=kにおける角速度センサのオフセット推定値gOKと、時刻T=k-1における角速度センサのオフセット推定値gOK-1とは等しい値であると仮定する。携帯機器1の角速度ωは、携帯機器1の利用者による携帯機器1の動かし方に依存して変化するため、時刻T=k-1の角速度ωk-1を用いて、時刻T=kにおける角速度ωkを定式化することは難しい。そこで、本実施形態に係る状態遷移モデルでは、便宜上、時刻T=kにおける角速度ωkと、時刻T=k-1における角速度ωk-1とは等しいと仮定する。前述の通り、磁気センサのオフセットqOFFは、携帯機器1の部品が発する内部磁界Biの方向及び大きさをセンサ座標系ΣSにおいて表現したベクトルである。従って、携帯機器1の内部状態が一定である間は、磁気センサのオフセットqOFFも変化しない。一方、携帯機器1の内部状態が変化した場合、例えば、携帯機器1に搭載された部品を流れる電流の大きさが変化した場合や、携帯機器1に搭載された部品の着磁状況が変化した場合には、磁気センサのオフセットqOFFも変化する。すなわち、携帯機器1の内部状態は、携帯機器1の利用者による携帯機器1の操作や、携帯機器1の外部の環境等に依存して変化する。従って、磁気センサのオフセットqOFFが変化するタイミングや磁気センサのオフセットqOFFの変化量を予測することは困難であり、時刻T=k-1における磁気センサのオフセット推定値qOK-1を用いて、時刻T=kにおける磁気センサのオフセット推定値qOKを定式化することは難しい。そこで、本実施形態に係る状態遷移モデルでは、便宜上、時刻T=kにおける磁気センサのオフセット推定値qOKと、時刻T=k-1における磁気センサのオフセット推定値qOK-1とは等しいと仮定する。このように、本実施形態に係る状態遷移モデルでは、以下に示す式27のように、状態ベクトルxkを構成する複数の状態変数のうち、姿勢μを表す状態変数以外は、前の時刻から変化しないという前提を置く。なお、カルマンフィルタの演算において、状態ベクトルxkは、観測残差zkに基づいて更新される。従って、式27に示された姿勢μ以外の値が、全く変化しないわけではない。状態ベクトルxkは、加算器129によって、観測残差zkに基づいて真値に近づくように更新される。具体的には、時刻T=kにおける状態ベクトルxkは、式6、式8、及び式10に示すように、時刻T=k-1における状態ベクトルxk-1を用いて算出される観測残差zkに基づいて更新される。同様に、時刻T=k-1における状態ベクトルxk-1は、時刻T=k-2における状態ベクトルxk-2を用いて算出される観測残差zk-1に基づいて更新される。すなわち、カルマンフィルタの演算において、状態ベクトルxkは、初期状態時刻T=0から時刻T=k-1に至る全ての状態ベクトルx0~xk-1の値を考慮して、更新される。従って、状態ベクトルxkの値が真値から大きく乖離した値に収束している場合、その後カルマンフィルタの演算を繰り返しても、状態ベクトルxkの値は真値に近づかなくなることがあり、また、仮に状態ベクトルxkの値が真値に近づくとしても、真値に近い値に収束するまでに長い時間を要することになる。前述の通り、磁気センサのオフセットqOFFは、携帯機器1の内部状態の変化に伴い、急激に大きく変化をする場合がある。磁気センサのオフセットqOFFが急激に大きく変化する場合、磁気センサのオフセット推定値qOKと、磁気センサのオフセットqOFF真値とは、大きく異なる値となる。磁気センサのオフセット推定値qOKが真値から大きく乖離した場合、カルマンフィルタによって状態ベクトルxkの更新を繰り返しても、磁気センサのオフセット推定値qOKが真値に近づかなくなることがあり、また、仮に真値に近づくことができても、真値に近い値に収束するまでに長い期間を要する。そして、磁気センサのオフセット推定値qOKが、真値とは大きく異なる値となる場合、携帯機器1は、地磁気Bgの向きや、携帯機器1の姿勢μ等を算出することができない。そこで、本実施形態では、中心点導出部200から中心点cOの座標が出力された場合、磁気センサのオフセット推定値qOを、中心点cOの座標によって上書きするステップS103の状態変数上書処理。中心点cOの座標は、N個の磁気データq1~qNに基づいて算出される。従って、中心点導出部200は、磁気センサのオフセットqOFFが急激に変化した場合であっても、オフセットqOFFの変化後に取得したN個の磁気データq1~qNを用いることで、磁気センサのオフセットqOFFを正確に表す中心点cOの座標を算出することができる。つまり、中心点cOの座標は、初期状態時刻T=0から時刻T=k-1に至る全ての状態ベクトルx0~xk-1の値を考慮して算出される磁気センサのオフセット推定値qOKに比べて、磁気センサのオフセットqOFFの急激な変化に柔軟に対応した真値に近い値として算出される。従って、本実施形態に係る状態推定装置は、磁気センサのオフセットqOFFが急激に変化した場合であっても、磁気センサのオフセット推定値qOを、中心点cOの座標によって上書きすることにより、状態ベクトルxkが真値から乖離することを防止することができる。[3.4.観測モデル]次に、観測モデル部125において行われる演算で用いる観測モデルについて説明する。3次元角速度センサ90から出力される角速度データgの推定値γgyroは、角速度ωと、角速度センサのオフセット推定値gOとを用いて、式28で与えられる。また、観測モデルにおいて、地上座標系ΣGでの地磁気ベクトルは式29で与えられる。従って、3次元磁気センサ70から出力される磁気データqの推定値γmagは、磁気センサのオフセット推定値qOと、式5で示した行列Bμとを用いて、式30で与えられる。また、観測モデルにおいて、地上座標系ΣGでの重力ベクトルは、重力加速度で正規化して、式31で与えられる。従って、3次元加速度センサ80から出力される加速度データaの推定値γaccは、式5で示した行列Bμを用いて、式32で与えられる。従って、式3、式28、式30、及び式32を、式8に適用することにより、観測残差zkは、式33で与えられる。[4.中心点導出部の動作]図7は、中心点導出部200の動作を表すフローチャートである。当該フローチャートに示される中心点導出処理は、CPU10が状態推定プログラムを実行することにより開始される。ステップS201において、CPU10は、初期化処理を実行する。具体的には、CPU10は、RAM20に記憶した磁気データを廃棄する。なお、本実施形態に係る初期化処理では、CPU10は、磁気データの全部を廃棄するが、本発明はこのような形態に限定するものではない。例えば、CPU10は、RAM20に蓄積されたデータのうち古い方から一定割合のデータのみを廃棄してもよい。ステップS202において、CPU10は、磁気データ蓄積処理を実行する。具体的には、CPU10は、3次元磁気センサ70から順次出力される磁気データqiをRAM20に格納する。そして、CPU10は、磁気データqiがRAM20にN個以上蓄積された場合、処理をステップS203へと進める。ステップS203において、CPU10は、磁気データ分布判定処理を実行する。具体的には、CPU10は、ステップS202の磁気データ蓄積処理においてRAM20等に蓄積された複数の磁気データのうち、新しい方からN個の磁気データq1~qNNは精度よく中心点cOの座標を導出するために必要な磁気データの規定測定回数を表す5以上の自然数が、後続するステップ204で実施される中心点算出処理において使用するデータとして適切なものであるか否かを判定する。そして、判定結果が肯定である場合、CPU10は、処理をステップS204へと進める。一方、判定結果が否定である場合、CPU10は、RAM20に蓄積された複数の磁気データの全部を廃棄した後に、処理をステップS202へと戻す。なお、本実施形態では、磁気データ分布判定処理における判定結果が否定である場合、CPU10は、RAM20に蓄積された複数の磁気データの全部を廃棄するが、本発明はこのような形態に限定されるものではない。例えば、CPU10は、判定結果が否定である場合に、RAM20に蓄積された複数の磁気データを廃棄せずに一定期間待機し、判定時にRAM20に蓄積されていた複数の磁気データに対して、待機している一定期間に新たにRAM20に蓄積された複数の磁気データを加えたうえで、再度判定処理を行ってもよい。ステップS204において、CPU10は、中心点算出処理を実行する。具体的には、CPU10は、ステップS202の磁気データ蓄積処理においてRAM20等に蓄積されたN個の磁気データq1~qNの示す座標が、ある球面第1球面Sの近傍に確率的に分布すると仮定して、球面Sの中心点cOの座標を算出する。ステップS205において、CPU10は、歪判定処理を実行する。具体的には、CPU10は、ステップS202の磁気データ蓄積処理においてRAM20等に蓄積されたN個の磁気データq1~qNの示す座標がある立体SDの表面近傍に分布すると仮定し、立体SDの形状と球面の形状とが同一と看做すことができる程度に近似しているか否かを判定する。判定結果が肯定の場合、すなわち、N個の磁気データq1~qNの示す座標が球面に近い形状の立体SDの表面近傍に分布すると判定した場合には、CPU10は、処理をステップS206に進める。一方、判定結果が否定の場合、すなわち、立体SDが球面とは異なる大きく歪んだ形状を有すると判定した場合、CPU10は、処理をステップS201に戻す。なお、本実施形態に係る歪判定処理では、立体SDと球面とが異なる形状であると判定した場合、CPU10は、処理をステップS201に戻すが、本発明はこれに限定されるものではない。例えば、歪判定処理において、立体SDと球面とが異なる形状であるとの判定結果が得られた場合に、CPU10は、何らかのメッセージを表示部50に出力した上で処理を一旦停止させ、携帯機器1の利用者からの指示を待ってステップS201から処理を再開させてもよい。N個の磁気データq1~qNを取得する際に、携帯機器1の利用者が携帯機器1を手で握り、携帯機器1の位置PSが変化するように携帯機器1を振るのではなく、携帯機器1の位置厳密には、3次元磁気センサ70の位置PSを固定したままその姿勢μのみを変化させるようにすると、外部磁界Bxの影響を低く抑えることができる図13B参照。そのため、ステップS205において、立体SDが球面とは異なる大きな歪を有する形状であると判定された場合には、携帯機器1の位置を固定したまま回転させることを、携帯機器1の利用者に対して指示してもよい。携帯機器1の利用者に対する指示は、携帯機器1の表示部50に画像や動画等を表示したり、音声等を用いたりすることで、行ってもよい。ステップS206において、CPU10は、中心点出力処理を実行する。具体的には、CPU10は、ステップS204の中心点算出処理で算出した中心点cOの座標を、前述した状態推定処理図5参照に引き渡す。CPU10は、以上に示したステップS201~S206を実行することにより、中心点導出部200として機能する。次にCPU10は、状態推定プログラムを終了する終了条件を充足するか否かを判定するステップS207。終了条件を充足する場合、CPU10は、当該フローチャートに示される中心点導出処理を終了する。一方、終了条件を充足しない場合、CPU10は、処理をステップS201に進める。なお、本実施形態では、CPU10は、ステップS207において、終了条件を充足しない場合に、処理をステップS201に進めるが、本発明はこのような形態に限定されるものではない。例えば、CPU10は、ステップS207の終了条件を充足しない場合、処理を一定期間停止させ、一定期間経過後に処理をステップS201に進めてもよい。また、CPU10は、ステップS207の終了条件を充足しない場合、状態ベクトル更新処理ステップS104で生成される観測残差zk等に基づいて中心点cOの座標の生成の要否を判定例えば、観測残差zkのノルムが閾値以上であるか否かを判定し、中心点cOの座標の生成が必要と判定された場合に、処理をステップS201へと進めてもよい。以下において、磁気データ分布判定処理ステップS203、中心点算出処理ステップS204、及び歪判定処理ステップS205の詳細について説明する。なお、説明の便宜上、まずステップS204の中心点算出処理について説明した後に、ステップS203の磁気データ分布判定処理について説明する。[5.中心点算出処理]中心点算出処理は、内部磁界Bi及び地磁気Bgの性質の違いを利用して、N個の磁気データq1~qNの示す座標から内部磁界Biの成分を分離することで、磁気センサのオフセットqOFFの候補となる座標である中心点cOの座標を算出する処理である。以下では、5.1節において、中心点算出処理の前提となる内部磁界Bi及び地磁気Bgの性質の相違点について説明したうえで、5.2節において、中心点算出処理の詳細について述べる。[5.1.内部磁界及び地磁気の性質]図8は、3次元磁気センサ70が、姿勢μをμ1~μNと変化させつつ内部磁界Bi及び地磁気Bgを検出した場合に3次元磁気センサ70から出力されるN個の磁気データq1~qNを、センサ座標系ΣSにおいてをプロットした図である。なお、図8に記載された各ベクトルの左上に付された添字Sは、当該ベクトルがセンサ座標系ΣSにおいて表されたベクトルであることを意味する。内部磁界Biは、携帯機器1の部品が発する磁界である。内部磁界Biは、携帯機器1の内部状態が変化しなければ、携帯機器1に対して一定の方向を向き、一定の大きさを有する磁界である。つまり、内部磁界Biは、センサ座標系ΣSにおいて、3次元磁気センサ70の位置Ps及び姿勢μに依存しない、一定の方向及び大きさを持つベクトルSBiとして表される。具体的には、図8に示すように、内部磁界を表すベクトルSBiは、センサ座標系ΣSの原点を起点とし、中心点cOGを終点とするベクトルである。一方、地磁気Bgは、磁極北に向かう水平成分と伏角方向の鉛直成分とを有する磁界であり、地上座標系ΣGにおいて一様な向き及び大きさを有する磁界である。よって、携帯機器1の姿勢μを変化させる場合は、地磁気Bgの大きさは変化しないが、携帯機器1から見た地磁気Bgの向きは変化する。つまり、地磁気Bgは、センサ座標系ΣSにおいて、携帯機器1の姿勢μに伴い方向を変化させる一定の大きさのベクトルSBgμとして表現される。複数の磁気データq1~qNの各々は、センサ座標系ΣSにおいて、内部磁界を表すベクトルSBiと、地磁気を表すベクトルSBgμとの和を表すベクトルの示す座標として出力される。従って、複数の磁気データq1~qNの各々がセンサ座標系ΣSにおいて示す座標は、中心点cOGを中心とし、地磁気Bgの大きさを半径とする球面SG上に分布する。なお、3次元磁気センサ70は、測定誤差を有するため、厳密には、複数の磁気データq1~qNの示す座標は球面SGの近傍に確率的に分布する。内部磁界を表すベクトルSBiが示す球面SGの中心点cOGの座標と、磁気センサのオフセットqOFFの示す座標とは、等しい。従って、磁気センサのオフセットqOFFは、複数の磁気データq1~qNに基づいて球面SGの中心点cOGの座標を求めることで算出される。つまり、磁気データqiの示す座標から、磁気センサのオフセットqOFFの示す座標を減算することにより、地磁気Bgの正確な方向及び大きさを算出することができる。[5.2.中心点算出処理の詳細]以下では、図9を参照しながら、ステップS204の中心点算出処理について説明する。中心点算出処理は、3次元磁気センサ70が出力するN個の磁気データq1~qNの示す座標が、半径rSの球面Sの近傍に分布すると仮定して、球面Sの中心点cOの座標を算出する処理である。なお、球面Sは、センサ座標系ΣS上で、複数の磁気データq1~qNの示す座標を近傍に有するように定められる球面の中心点を求めるために、計算の便宜上導入された球面であり、地磁気Bgを表す球面SGとは異なる球面である。しかし、3次元磁気センサ70が、内部磁界Bi及び地磁気Bgのみを検出する場合には、複数の磁気データq1~qNの示す座標は球面SGの近傍に分布することになる。従って、この場合、複数の磁気データq1~qNの示す座標を近傍に有するように定められる球面Sと球面SGとは、ほぼ同一の球面となり、球面Sの中心点cOと球面SGの中心点cOGとは、ほぼ同一の座標となる。よって、3次元磁気センサ70が、内部磁界Bi及び地磁気Bgのみを検出する場合、複数の磁気データq1~qNの示す座標を近傍に有するように定められる球面Sの中心点cOの座標を算出することで、球面SGの中心点cOGの座標、すなわち、磁気センサのオフセットqOFFを算出することができる。以下、具体的に、球面Sの中心点cOの座標の算出方法について述べる。なお、以下において記載されるベクトル及び座標は、特に断りが無い場合には、センサ座標系ΣSにおいて表現されたものであるとする。複数の磁気データq1~qNの示す座標が、半径rSの球面S上に存在すると仮定する場合、各磁気データqiの示す座標と、球面Sの中心点cOとの距離はrSであるため、以下の式34が成立する。なお、磁気データqiの示す座標の各要素は式35で表され、中心点cOの座標の各要素は式36で表される。ここで、複数の磁気データq1~qNの示す座標の重心qCを、式37で定義する。なお、重心qCの各要素は、式38で示される。磁気データqiの示す座標及び中心点cOの座標を、重心qCを起点とするベクトルとして表現することで、式34は、式39に変形される。式39の磁気データqiの示す座標に対して、複数の磁気データq1~qNの示す座標の各々を代入した結果を、磁気データの個数Nで割り算することにより、以下の式40が得られる。そして、式39と式40との差分を取ることで、以下の式41が得られる。このように、磁気データである変数qiの示す座標が、球面S上に存在することを示す式34は、式41に変形された。この式41の変数qiに対して、N個の磁気データq1~qNの示す座標をそれぞれ代入することで得られるN個の方程式は、行列Xを用いて式42として表現される。なお、行列Xは、式43に示すN行3列の行列である。また、ベクトルjは、式44に示すN次元のベクトルであり、値Rは、式45に示す値である。式42は、複数の磁気データq1~qNの示す座標の全てが、中心点cOを中心とする球面S上に完全に一致する場合には解を有する。しかし、3次元磁気センサ70の測定誤差等を考慮すると、複数の磁気データq1~qNの示す座標の全てが、球面Sと完全に一致することは無いため、式42は解を持たない。そこで、統計的な手法により尤もらしい解を得るために、式46で表す誤差を吸収するベクトルである第1球面誤差ベクトルδSを導入する。ここで、式47により表される3次元のベクトルcは、中心点cOの座標を表すための変数ベクトルである。第1球面誤差ベクトルδSのノルムを最小にするベクトルc、換言すれば、Tを最小化するようなベクトルcにより示される座標が、球面Sの中心点cOの座標として尤もらしいものであるいえる。ここで、以下の式48で示される目的関数中心点算出関数fScを定義すると、目的関数fScを最小化するベクトルcにより示される座標が、球面Sの中心点cOの座標として尤もらしい値となる。中心点cOの座標は、式50に示す共分散行列Aが正則である場合には、式49で表される。前述の通り、3次元磁気センサ70が、内部磁界Bi及び地磁気Bgのみを検出する場合には、球面Sと、地磁気Bgを表す球面SGとは、ほぼ同一の球面となり、球面Sの中心点cOと、球面SGの中心点cOGとは、ほぼ同一の座標となる。従って、3次元磁気センサ70が内部磁界Bi及び地磁気Bgのみを検出する場合、式49に示される中心点cOの座標を示すベクトルを、磁気センサのオフセットqOFFとして採用することができる。[6.磁気データ分布判定処理]以下では、ステップS203で行われる、磁気データ分布判定処理について説明する。前述した、中心点算出処理において、球面Sの中心点cOを算出するためには、複数の磁気データq1~qNの示す座標が、センサ座標系ΣSにおいて3次元的な広がりを有して分布していることが必要となる。しかし、携帯機器13次元磁気センサ70の姿勢μは、携帯機器1の利用者が手で携帯機器1を動かすことにより変化するため、携帯機器1の動かし方が不十分な場合には、携帯機器1の姿勢変化は3次元的とはならず2次元的なものとなることがある。この場合、センサ座標系ΣSにおける複数の磁気データq1~qNの示す座標は、3次元的な広がりを有さず2次元的に分布する。例えば、図10に示すように、複数の磁気データq1~qNの示す座標が、センサ座標系ΣSの平面π上の円πC近傍に2次元的に分布する場合、球面Sは、円πCを切断面に有するような球面であるということしか特定できない。円πCを切断面に有するような球面は、円πCの中心点πCOを通り平面πに直交する直線πL上の中心点cπ1を中心とする球面Sπ1であるかもしれないし、直線πL上の中心点cπ2を中心とする球面Sπ2であるかもしれない。つまり、球面Sの中心点cOは、直線πL上に位置することまでは特定可能であるが、具体的に直線πL上のどの位置に存在するかについては特定不可能である。このように、複数の磁気データq1~qNの示す座標が2次元的に分布する場合、複数の磁気データq1~qNに基づいて、正確な中心点cOの座標を算出することはできない。複数の磁気データq1~qNに基づいて球面Sの中心点cOの座標を算出するためには、図11に示すように、センサ座標系ΣSにおいて複数の磁気データq1~qNの示す座標が3次元的に広がりを持って分布していることが必要となる。磁気データ分布判定処理において、CPU10中心点導出部200は、複数の磁気データq1~qNの示す座標が3次元的に分布しているか否かを、式50に示される共分散行列Aを用いて判定する。以下において、共分散行列Aの性質を説明する。共分散行列Aの固有値を、大きい順番に、最大固有値λ1、中間固有値λ2、最小固有値λ3とし、それぞれの固有値に対応する大きさ1に正規化された固有ベクトルをu1、u2、u3とする。また、磁気データqiの示す座標を、前述した重心qCを原点とする重心座標系ΣCにおいて表したベクトルをCqiと表す。このとき、固有値λjj=1、2、3は、固有ベクトルuj方向の分散σ2jに等しい。図11に示すように、各固有ベクトルu1~u3を重心座標系ΣCの原点qCを起点となるように配置する。このとき、例えばj=1の場合について検討する。固有値λ1は、ベクトルCqiを、固有ベクトルu1へ射影した長さLi1の二乗Li12を、N個の磁気データCqii=1、2、...Nについて平均した値に等しい。すなわち、固有値λjは、N個の磁気データqiの示す座標が、重心qCから固有ベクトルuj方向にどの程度離れているか、つまり、複数の磁気データq1~qNの示す座標の分布が固有ベクトルujの方向にどの程度の広がりを有しているかを表す値である。最小固有値λ3に対応する固有ベクトルu3の方向が、複数の磁気データq1~qNの示す座標の分布の広がりが最も小さい方向であり、最小固有値λ3が、複数の磁気データq1~qNの示す座標の分布の広がりが最も小さい方向における広がりの程度を示す指標である。従って、複数の磁気データq1~qNの示す座標が3次元的に分布しているといえるためには、最小固有値λ3の値が、一定の閾値分散許容値λ0以上であればよい。磁気データ分布判定処理において、中心点導出部200は、共分散行列Aの最小固有値λ3が閾値λ0以上であれば、複数の磁気データq1~qNの示す座標の分布が十分に3次元的であると判断して、処理を前述したステップS204の中心点算出処理へと進める。一方、中心点導出部200は、最小固有値λ3が閾値λ0未満である場合には、複数の磁気データq1~qNの示す座標が3次元的な広がりを有さず2次元的に分布すると判断し、処理をステップS202の磁気データ蓄積処理に戻す。[7.歪判定処理]ところで、3次元磁気センサ70は、図3に示したように、地磁気Bg及び内部磁界Biの他に、携帯機器1の外部に存在する物体2が出す外部磁界Bxを検出する場合がある。外部磁界Bxは、物体2と携帯機器1との相対的位置関係により、方向及び大きさが変化する不均一な磁界である。5節で説明した中心点算出処理は、3次元磁気センサ70が、地磁気Bg及び内部磁界Biのみを検出することを前提として、磁気センサのオフセットqOFFの示す座標の候補である球面Sの中心点cOの座標を算出するものである。従って、3次元磁気センサ70が、地磁気Bg及び内部磁界Biの他に、不均一な外部磁界Bxを検出する場合、中心点算出処理により算出される中心点cOの座標は、外部磁界Bxの影響による誤差を有し、内部磁界Biすなわち、磁気センサのオフセットqOFFを表す座標とは異なる座標となる可能性が高い。このような中心点cOの座標を示すベクトルを磁気センサのオフセットqOFFとして採用した場合、正しい補正処理を行うことができず、携帯機器1は正確な地磁気Bgの方向及び大きさを算出することができない。ステップS205の歪判定処理において、CPU10中心点導出部200は、不均一な外部磁界Bxの有無、及び外部磁界Bxの影響の大きさを評価することにより、中心点算出処理により算出された中心点cOの座標を示すベクトルを、磁気センサのオフセットqOFFとして採用することの可否を判断する。以下では、7.1節において、歪判定処理における判定の対象となる外部磁界Bxの性質について説明したうえで、7.2節において、歪判定処理の詳細を説明する。[7.1.外部磁界の性質]図12は、3次元磁気センサ70の位置PsをPS1~PSNと変化させると共に、姿勢μをμ1~μNと変化させて磁界を測定したときの、3次元磁気センサ70が出力する複数の磁気データq1~qNの示す座標を、センサ座標系ΣSにおいてプロットした図である。なお、図12では、内部磁界Bi、地磁気Bg、及び外部磁界Bxが存在する場合を想定している。外部磁界Bxは不均一な磁界であり、外部磁界Bxを発生させる物体2と、携帯機器1との相対的位置関係により、方向及び大きさが変化する。従って、3次元磁気センサ70の位置Psが変化した場合、3次元磁気センサ70が検出する外部磁界Bxの向き及び大きさも変化する。同様に、3次元磁気センサ70の姿勢μが変化した場合、3次元磁気センサ70が検出する外部磁界Bxの向きも変化する。すなわち、外部磁界Bxは、センサ座標系ΣSにおいて、3次元磁気センサ70の姿勢μ及び位置Psに依存して方向及び大きさの双方を変化させるベクトルSBxμ、Psとして表現される。複数の磁気データq1~qNの示す座標の各々は、内部磁界を表すベクトルSBi、地磁気を表すベクトルSBgμ、及び外部磁界を表すベクトルSBxμ、Psの和を表すベクトルにより示される。従って、複数の磁気データq1~qNの示す座標は、中心点cOGを起点とする地磁気を表すベクトルSBgμの終点を表す球面SGと、中心点cOGを起点とする外部磁界を表すベクトルSBxμ、Psの終点を表す曲面SXとを、重ね合わせた立体SDの表面近傍に分布する。外部磁界Bxを表す曲面SXが、球面とは異なる歪んだ形状を有する場合には、立体SDも球面とは異なる歪んだ形状を有する。複数の磁気データq1~qNの示す座標が、球面とは異なる歪んだ形状の立体SDの表面近傍に分布する場合、複数の磁気データq1~qNの示す座標をできるだけ近傍に有するように球面Sを定めても、複数の磁気データq1~qNの示す座標と球面Sとの間の距離は大きく離れる。すなわち、この場合は、立体SDの表面近傍に分布する複数の磁気データq1~qNの示す座標は、球面Sの近傍には分布しない。同様に、磁気データq1~qNの示す座標は、地磁気Bgを表す球面SGの近傍にも分布しない。従って、複数の磁気データq1~qNの示す座標が、球面とは異なる歪んだ形状の立体SDの表面近傍に分布する場合、球面Sの中心点cOと球面SGの中心点cOGとは異なる座標となる可能性が高く、中心点cOの座標を示すベクトルを、磁気センサのオフセットqOFFとして採用することはできない図15参照。但し、不均一な外部磁界Bxの影響が小さく、立体SDの形状がほぼ球面と看做せる場合、複数の磁気データq1~qNの示す座標を近傍に有するように定められる球面Sと、立体SDとをほぼ同一の図形と看做すことが可能となり、また、球面Sの中心点cOと、地磁気を表す球面SGの中心点cOGとはほぼ同一の座標となる。従って、立体SDの形状がほぼ球面と看做せる場合は、球面Sの中心点cOの座標を示すベクトルを、磁気センサのオフセットqOFFとして採用することができる。以下、図13を参照しながら、外部磁界Bxの影響が小さい場合について具体的に説明する。図13Aは、外部磁界Bxが微弱である場合を示す。立体SDは、地磁気Bgを表す球面SGと、外部磁界Bxを表す曲面SXとの重ね合わせである。従って、外部磁界Bxが微弱であり、外部磁界Bxを表す曲面SXが小さい場合には、立体SDと球面SGとは、ほぼ同一の図形となる。このとき、立体SDの表面近傍に分布する複数の磁気データq1~qNの示す座標は、球面SGの近傍にも分布するものと看做すことができる。また、複数の磁気データq1~qNの示す座標を近傍に有する立体SDの形状は、ほぼ球面であるため、複数の磁気データq1~qNの示す座標を近傍に有するように定められる球面Sと、立体SDとは、ほぼ同一の図形となる。すなわち、立体SD、球面SG、及び球面Sは、ほぼ同一の図形になり、球面Sの中心点cOと球面SGの中心点cOGとは、ほぼ等しい座標になる。従って、図13Aに示すような、外部磁界Bxが微弱である場合には、中心点算出処理により算出された中心点cOの座標を示すベクトルを、磁気センサのオフセットqOFFとして採用することができる。なお、図13Bに示すように、不均一な外部磁界Bxが大きい場合であっても、立体SDの形状がほぼ球面と看做せる場合がある。例えば、不均一な外部磁界Bxが存在する場合であっても、N個の磁気データq1~qNを取得する際に、携帯機器1の利用者が携帯機器1を手で握って携帯機器1の位置PSが変化するように振るのではなく、携帯機器1の位置PS厳密には、3次元磁気センサ70の位置を固定したままその姿勢μのみを変化させる場合には、外部磁界Bxは、センサ座標系ΣSにおいて、姿勢μに基づいてその方向のみを変化させる一定の大きさのベクトルSBxμとして表現される。この場合、外部磁界Bxを表す曲面SXの形状は、中心点cOGを中心とする球面となるため、球面SGと曲面SXとの重ね合わせである立体SDの形状もほぼ球面となる。また、立体SDは、中心点cOGを中心とする球面SGと、中心点cOGを中心とする球面の形状を有する曲面SXとを重ね合わせた図形であるため、立体SDが表す球面の中心点は、球面SGの中心点cOGと、ほぼ一致する。前述の通り、立体SDの形状がほぼ球面である場合には、複数の磁気データq1~qNの示す座標を近傍に有するように定められる球面Sと、立体SDとは、ほぼ同一の図形となる。よって、球面Sの中心点cOと、球面SGの中心点cOGとは、ほぼ等しい座標になる。従って、図13Bに示すような、外部磁界Bxが一定の大きさのベクトルとして検出される場合には、中心点算出処理により算出された中心点cOの座標を示すベクトルを、磁気センサのオフセットqOFFとして採用することができる。このように、不均一な外部磁界Bxの影響が大きい場合、すなわち、立体SDが球面とは異なる歪んだ形状を有する場合には、中心点算出処理により算出された中心点cOの座標を示すベクトルを、磁気センサのオフセットqOFFとして採用することはできない。一方、不均一な外部磁界Bxが微弱な場合、または、3次元磁気センサ70が不均一な外部磁界Bxを一定の大きさのベクトルとして検出する場合等、立体SDの形状がほぼ球面と看做すことができる場合には、球面Sの中心点cOの座標を示すベクトルを磁気センサのオフセットqOFFとして採用することができる。[7.2.歪判定処理の詳細]以下では、図13乃至図15を参照しながら、ステップS205の歪判定処理の詳細について説明する。前述した中心点算出処理では、複数の磁気データq1~qNの示す座標が、式34で示される球面Sの近傍に分布すると仮定した。これに対して本節で説明する歪判定処理は、複数の磁気データq1~qNの示す座標が、球面Sとは異なる歪んだ形状の立体SDの表面近傍に分布すると仮定する。立体SDは、式51で表される。式51に示される立体SDは、図14に示すように、球面を表す成分Xc-qC-jと、式52に示すN次元ベクトルの歪誤差ベクトルkEとを、加算したものである。なお、式51の右辺に現れる0Nの右下の添え字Nは、零ベクトルの次元を表している。また、式53に示す歪評価行列Eは、3行3列の対称行列である。なお、式51に現れる成分Xc-qC-jが表す球面と、中心点算出手段により算出される球面Sとは、必ずしもセンサ座標系ΣS上で同一の図形を表すものではない。そこで、以下では、式51に現れる成分Xc-qC-jとして表される球面を球面第2球面S2と称し、球面S2の中心点を中心点第2球面の中心点cO2と称する。また、式51において用いられるベクトルcは、式47に示した3次元のベクトルであるが、本節では、中心点cO2の示す座標を算出するための変数ベクトルとして用いる。歪判定処理は、式51のうち、歪を表す成分kEの大きさを評価することにより、立体SDの形状と、球面S2の形状とが、どの程度相違しているかについて評価するものである。歪誤差ベクトルkEを構成する要素のうち、第i行目の要素keは、以下の式54で表される関数kepに対して、中心点cOを起点として磁気データqiを表すベクトルqi-cOを代入することにより与えられる。関数kepは、式53に示す歪評価行列Eを係数行列とし、式55に示すベクトルpの3つの要素を変数とする、二次形式で表現される関数である。つまり、関数kepは、ベクトルpと、ベクトルpを歪評価行列Eにより変換したベクトルEpとの内積を示す。従って、歪誤差ベクトルkEのうち、第i行目の要素keqi-cOは、球面Sの中心点cOを起点として磁気データqiの座標を表したベクトルqi-cOと、ベクトルqi-cOを歪評価行列Eにより変換して得られるベクトルEqi-cOとの内積である。なお、3次元磁気センサ70の測定誤差等を考慮すると、複数の磁気データq1~qNの各々が示す座標の全てが、立体SDと完全に一致する位置に存在することは無いため、式51は解を持たない。そこで、統計的な手法により尤もらしい解を得るために、式56で表される誤差を吸収するベクトルである立体誤差ベクトルδSDを導入する。立体誤差ベクトルδSDは、第2球面誤差ベクトルδS2と、歪誤差ベクトルkEとを、加算したものである。なお、詳細は後述するが、第2球面誤差ベクトルδS2は、式51のうち、球面を表す成分Xc-qC-jに対応する成分である。立体誤差ベクトルδSDは、複数の磁気データq1~qNの座標と、立体SDの表面との誤差を表す、N次元のベクトルである。立体誤差ベクトルδSDのノルムを最小にするようなベクトルc及び歪評価行列E、換言すれば、57で示される目的関数歪評価関数fSDE、cを最小化するようなベクトルc及び歪評価行列Eによって、複数の磁気データq1~qNの示す座標を表面近傍に有する立体SDが表現される。歪判定処理は、立体SDの形状における歪誤差ベクトルkEの影響の大きさを、後述する式58及び59で示される歪評価値gDEに基づいて評価する処理である。立体SDの形状における歪誤差ベクトルkEの影響の大きさを評価することで、立体SDと球面Sとの形状の相違の程度について把握することが可能となる。以下で、式56に示した立体誤差ベクトルδSDの性質を、式46に示した第1球面誤差ベクトルδSの性質と対比しつつ説明する。まず、第1球面誤差ベクトルδSは、複数の磁気データq1~qNが示す座標と、球面Sとの誤差を吸収するためのベクトルである。第1球面誤差ベクトルδSを構成する1行目~N行目の各要素は、それぞれ独立した変数である。従って、第1球面誤差ベクトルδSによって、複数の磁気データq1~qNが示す座標と球面Sとの誤差を吸収する場合には、複数の磁気データq1~qNが示す座標と球面SとのN個の誤差の各々は、互いに制約の無い、N個が互いに独立して定められる値となる。つまり、第1球面誤差ベクトルδSによって表されるN個の誤差は、それぞれが独立に確率的に定められるものであり、N個の誤差は全体として、対称性を有し、かつ方向依存性の無いホワイトノイズである。すなわち、中心点算出処理は、複数の磁気データq1~qNが示す座標と、球面Sとの誤差をホワイトノイズである第1球面誤差ベクトルδSにより表現し、第1球面誤差ベクトルδSを最小化するような球面Sの中心点cOの座標を求める処理である。一方、立体誤差ベクトルδSDは、第2球面誤差ベクトルδS2と歪誤差ベクトルkEとの和によって表されるベクトルであり、磁気データq1~qNが示す座標と立体SDとの誤差を吸収するベクトルである。第2球面誤差ベクトルδS2は、第1球面誤差ベクトルδSと同様に式46の右辺で表され、磁気データq1~qNの座標と、球面S2との誤差を、ホワイトノイズとして表現するベクトルである。一方、歪誤差ベクトルkEは、式54で示した3変数二次形式の関数kepを各要素とするベクトルである。3変数の二次形式は、変数が2次の項から構成される関数であり、3次元空間上の様々な曲面、例えば、直線、平面、柱面、球面、楕円面、錐面、1葉双曲面、2葉双曲面、及び各種放物面等を描くことができる。従って、歪誤差ベクトルkEは、各磁気データq1~qNが示す座標と球面S2とのN個の誤差の各々を、互いに独立したものではなく、N個の誤差の全てが同一の関数kepにより表される3次元空間上の曲面上に存在するという制約を持った値として表現する。このように、立体誤差ベクトルδSDは、各磁気データq1~qNが示す座標と球面S2とのN個の誤差を、ホワイトノイズである第2球面誤差ベクトルδS2と、滑らかな曲面による球面S2からの歪という性質を有する歪誤差ベクトルkEとに分離して表現する。以下では、立体SDの形状における歪誤差ベクトルkEの影響の大きさと、中心点算出処理で算出された球面Sの中心点cOの座標を示すベクトルを磁気センサのオフセットqOFFとして採用することの可否についての判断との関係について、図13及び図15を参照しつつ説明する。図13A及びBは、立体SDの形状における歪誤差ベクトルkEの影響が小さい場合を示している。歪誤差ベクトルkEの影響が小さい場合とは、歪誤差ベクトルkEの大きさが小さいため、立体SDを表す式51が成立すれば、同時に、球面Sを表す式42についても成立するものと看做せる場合を意味する。すなわち、歪誤差ベクトルkEの影響が小さい場合には、式57に定める目的関数fSDE、cを最小化して求められる立体SDと、式48に定める目的関数fScを最小化して求められる球面Sとが、ほぼ等しくなる。前述の通り、立体SDの形状が球面と看做せる場合には、複数の磁気データq1~qNが示す座標は、立体SDの表面近傍に分布すると同時に、球面Sの近傍にも分布すると看做すことができる。そして、立体SDが表す球面の中心点は、球面SGの中心点cOGと、ほぼ一致するため、球面Sの中心点cOと球面SGの中心点cOGとがほぼ等しい座標となる。従って、歪誤差ベクトルkEの影響が小さい場合には、球面Sの中心点cOの座標を示すベクトルを磁気センサのオフセットqOFFとして採用することができる。なお、歪評価行列Eが零行列となる場合には、歪誤差ベクトルkEも零ベクトルとなるため、式48と式57とが等しい式となり、立体SDと球面Sとは完全に同一の図形となる。図15は、歪誤差ベクトルkEの影響が大きい場合を示している。歪誤差ベクトルkEの影響が大きい場合とは、複数の磁気データq1~qNが示す座標をできるだけ近傍に有するように球面S2を定めても、複数の磁気データq1~qNが示す座標と球面S2との誤差をホワイトノイズである第2球面誤差ベクトルδS2によって吸収することができず、球面S2からの滑らかな歪を表す歪誤差ベクトルkEにより吸収する必要がある場合を示している。すなわち、この場合、目的関数fSDE、cを最小化して求められる立体SDは、目的関数fScを最小化して求める球面Sとは大きく異なる歪んだ形状になる。前述の通り、立体SDの形状が球面とは異なる歪んだ形状を有する場合には、立体SDの表面近傍に分布する複数の磁気データq1~qNが示す座標は、球面Sの近傍に分布すると看做すことはできない。中心点算出処理は、球面Sの近傍に磁気データq1~qNの示す座標が存在していることを前提として、磁気センサのオフセットqOFFの候補である中心点cOの座標を算出する処理であるため、磁気データq1~qNの示す座標を近傍に有さない球面Sの中心点cOは、磁気センサのオフセットqOFFとしての意味を有さない。従って、この場合、中心点cOの座標を示すベクトルを磁気センサのオフセットqOFFとして採用することはできない。このように、立体SDの形状における歪誤差ベクトルkEの影響が小さい場合、立体SD上に分布する複数の磁気データq1~qNの示す座標は、球面S上にも分布すると看做すことができ、中心点算出手段で算出した球面Sの中心点cOの座標を示すベクトルを磁気センサのオフセットqOFFとして採用することができるが、立体SDの形状における歪誤差ベクトルkEの影響が大きければ、中心点cOの座標は誤差を有する不正確な値となり、中心点cOの座標を示すベクトルを磁気センサのオフセットqOFFとして採用することを防止する必要がある。以下では、立体SDの形状における歪誤差ベクトルkEの影響の大きさを評価する方法について述べる。ここで、以下の式58及び59に示す歪評価値gDEを、立体SDの形状における歪誤差ベクトルkEの影響の大きさを評価する評価値として定義する。歪評価値gDEは、歪評価行列Eの有する3つの固有値のうち、絶対値が最大となる最大固有値λE1の絶対値歪評価行列Eのノルムである。歪評価値gDEの値が一定の閾値歪許容値δ0以下の小さな値であれば、立体SDと球面S2とは等しい図形であると看做すことができるため、立体SDの表面近傍に分布する複数の磁気データq1~qNの示す座標は、球面Sの近傍に分布すると看做すこともでき、中心点算出処理により求められた球面Sの中心点cOの座標を示すベクトルを磁気センサのオフセットqOFFとして採用することができる。以下に、歪評価値gDEを求める方法について説明する。まず、式54に示した関数kepは、二次形式で表される関数であるため、以下の式60のように、ベクトルの内積として表現することもできる。また、歪誤差ベクトルkEの第i行目の要素keは、関数kepに対して、ベクトルを代入した値であるため、式62によって示される6次元のベクトルke2i、及び、式63によって表される歪評価行列Eの各成分を並べた6次元のベクトルeEを用いて、式61のように表せる。ここで、式64で示す行列X2を導入する。行列X2は、各行に、式62によって示される6次元のベクトルke2iを転置したものと、3次元のベクトルを転置したものを並べた、N行9列の行列である。式64で表される行列X2を用いることで、式57に示される目的関数fSDE、cを、式65で表される目的関数gSDeに変形することができる。なお、ベクトルeは、式66に示すように、式63に示した6次元のベクトルeEと、式67に示す3次元のベクトルeXとを並べた、9次元のベクトルである。また、式67に示すベクトルeXは、式47に示す変数であるベクトルcから、式37に示した重心qCを引いたベクトルである。式65で表される目的関数gSDeを最小化する解e=e0は、以下の式68式で表現される連立方程式に対して、ガウス消去法や、コレスキー分解法などを適用することで、求められる。なお、式68は、式65に対して、最小二乗法を適用することによって算出される連立方程式である。このようにして得た解e0に基づき、式53の歪評価行列Eを復元する。そして、式58に示した歪評価値gDEの値、即ち、歪評価行列Eのノルムを求め、歪評価値gDEの値が閾値δ0以下であるか否かを判定する。なお、歪評価行列Eのノルムは、歪評価行列Eの有する3つの固有値のうち絶対値が最大となる固有値λE1の絶対値に等しいため、ヤコビ法や、冪乗法により求めることができる。そして、歪評価値gDEの値が閾値δ0以下である場合、中心点導出部200は、処理をステップS206の中心点出力処理に進める。一方、歪評価値gDEの値が閾値δ0よりも大きな値である場合、中心点導出部200は、処理をステップS201の初期化処理に戻す。[8.中心点出力処理]ステップS206の中心点出力処理は、歪判定処理において歪評価値gDEの値が閾値δ0以下であると判定された場合に実行される処理である。中心点出力処理において、CPU10中心点導出部200は、中心点算出処理において算出された球面Sの中心点cOの座標を出力し、状態推定部100に対して供給する。本実施形態では、中心点導出部200は、球面Sの中心点cOの座標を出力しているが、球面S2の中心点cO2の座標を出力してもよい。歪評価値gDEの値が閾値δ0以下である場合には、中心点算出処理により求められた球面Sの中心点cOと、球面S2の中心点cO2とは、ほぼ等しい座標となり、中心点cOの座標を示すベクトルまたは中心点cO2の座標を示すベクトルのいずれも、磁気センサのオフセットqOFFとして採用することができるからである。なお、球面S2の中心点cO2の座標は、目的関数gSDeを最小化する解e0のうち式66のeXに相当する部分を式67に代入することにより求められる。[9.結論]以上に示したように、本実施形態に係る状態推定装置は、中心点導出部200と、カルマンフィルタ演算部120を備える状態推定部100とを備え、カルマンフィルタ演算部120がカルマンフィルタの演算に用いる状態ベクトルxkの要素のうち磁気センサのオフセット推定値qOを、中心点導出部200が出力する中心点cOの座標によって上書きをしてカルマンフィルタの演算を行う。前述の通り、カルマンフィルタの状態遷移モデルにおいて磁気センサのオフセット推定値qOの経時的な変化を定式化することはできず、また、カルマンフィルタによる状態ベクトルxkは更新前の状態ベクトルxkの影響を受けた値となるため、磁気センサのオフセット推定値qOは、真値磁気センサのオフセットqOFFから大きく乖離した値となる場合がある。この場合、状態推定装置は、システムの状態を正確に推定できないため、状態推定装置が算出する状態変数に基づいて、地磁気Bgの向きや携帯機器1の姿勢μ等を算出することはできない。これに対して、本実施形態におけるカルマンフィルタ演算部120は、中心点導出部200が出力した中心点cOの座標により、状態ベクトルxkの要素である磁気センサのオフセット推定値qOを上書きするステップS103において行われる、状態変数上書処理。前述の通り、中心点cOの座標は、比較的短い期間に3次元磁気センサ70から出力された磁気データq1~qNに基づいて算出される。すなわち、中心点cOの座標を示すベクトルは、磁気センサのオフセットqOFFが急激且つ大きな変化をした場合であっても、変化後の磁気センサのオフセットqOFFを正確に表すことができる。このような中心点cOの座標によって、磁気センサのオフセット推定値qOを上書きすることにより、状態ベクトルxkが真値から大きく異なる値となることを防止した。そして、携帯機器1が、状態ベクトルxkに基づいて、正確な地磁気Bgの方向及び大きさ等を算出することを可能とした。また、本実施形態に係る、中心点導出部200は、中心点算出処理を行う前に、磁気データ分布判定処理を行うことで、中心点算出処理で用いる複数の磁気データq1~qNが、中心点算出処理を行う上で適切なデータであるか否かの判定を行う。これにより、複数の磁気データq1~qNの示す座標の分布の3次元的な広がりの程度を把握し、複数の磁気データq1~qNの示す座標が平面的に分布している場合には、中心点算出処理を行わない。この結果、平面的に分布する複数の磁気データq1~qNに基づいた、不正確な真値とは異なる中心点cOの座標の算出を未然に防止することが可能となり、不要な演算処理の防止による低消費電力化が可能となる。また、本実施形態に係る、中心点導出部200は、中心点算出処理を行った後に、歪判定処理を行い、中心点算出処理で算出された球面Sの中心点cOの座標が、外部磁界Bxによる影響を受けた不適切な真値とは異なる値であるか否かを判定する。これにより、外部磁界Bxによる影響を受けた不適切な磁気センサのオフセットqOFFと看做すことのできない中心点cOの座標の出力を防止することが可能となる。このように、本実施形態に係る中心点導出部200は、磁気データ分布判定処理及び歪判定処理により、磁気センサのオフセットqOFFと看做すことのできる真値に近い中心点cOの座標のみを出力する。そして、カルマンフィルタ演算部120は、真値に近い中心点cOの座標により磁気センサのオフセット推定値qOを上書きすることで、状態ベクトルxkをより真値に近い正確な値へと更新することが可能となる。従って、本実施形態に係る携帯機器1は、状態推定装置が算出する状態ベクトルxkに基づいて、正確な地磁気Bgの値を安定的に算出することが可能となった。<B.変形例>本発明は上述した実施形態に限定されるものではなく、例えば、以下の変形が可能である。また、以下に示す変形例のうちの2以上の変形例を組み合わせることもできる。1変形例1上述した実施形態では、カルマンフィルタ演算部120は、状態ベクトルxのうち磁気センサのオフセット推定値qOを表す要素を、中心点導出部200から出力された中心点cOの座標によって上書きをしてカルマンフィルタの演算を行うものであったが、本発明はこのような形態に限定されるものでは無く、状態ベクトルxを中心点cOの座標によって上書きせずにカルマンフィルタの演算を実行してもよい。図16は、変形例1に係る状態推定装置のうち、CPU10が変形例1に係る状態推定プログラムを実行することにより実現される機能を表す機能ブロック図である。図16に示すように、変形例1に係る状態推定装置は、出力統合部300を備える点と、状態推定部100の代わりに状態推定部100aを備える点とを除いて、実施形態に係る状態推定装置と同様に構成される。状態推定部100aは、カルマンフィルタ演算部120の代わりにカルマンフィルタ演算部120aを備える点を除き、状態推定部100と同様に構成される。カルマンフィルタ演算部120aは、状態ベクトルxを中心点cOの座標によって上書きしない点を除き、カルマンフィルタ演算部120と同様に動作する。出力統合部300は、カルマンフィルタ演算部120aが出力する更新後の状態ベクトルxkと、中心点導出部200が出力する中心点cOの座標とに基づいて、状態ベクトルxckを出力する。状態ベクトルxckは、状態ベクトルxkと同様に、姿勢μ、地磁気の強さr、地磁気の伏角φ、角速度ω、角速度センサのオフセット推定値gO、及び磁気センサのオフセット推定値qOを要素とするベクトルである。状態ベクトルxckは、初期状態時刻k=0から、中心点導出部200が中心点cOの座標を出力するまでの期間は、状態ベクトルxkと等しい値に設定される。一方、中心点導出部200が中心点cOの座標を出力した後は、状態ベクトルxckを構成する要素のうち、磁気センサのオフセット推定値qOは、中心点cOの座標と等しい値に設定される。カルマンフィルタ演算部120aは、3次元地磁気センサ70、3次元加速度センサ80、及び3次元角速度センサ90の3つのセンサの出力を統合して、各種状態変数を推定するものであるため、高速に磁気センサのオフセット推定値qOを算出することができる。一方、中心点導出部200は、N個の磁気データq1~qNに基づいて中心点cOの座標を算出するため、中心点cOの座標が算出されるまでに一定の時間を要するが、歪判定処理により中心点cOの座標の良否判定を行っているため、中心点cOの座標を示すベクトルは磁気センサのオフセットqOFFと看做すことができる程度に正確な値である。従って、変形例1に係る携帯機器1は、中心点導出部200が中心点cOの座標を算出するまでの間は、カルマンフィルタ演算部120aが算出した磁気センサのオフセット推定値qOを用いて地磁気Bgのおおよその方向及び大きさを得ることができる。これにより、携帯機器1は、携帯機器1の利用者を長時間待たせること無く、地磁気Bgの方向及び大きさを表示することが可能となり、利便性が向上する。2変形例2上述した実施形態では、中心点導出部200は、初期化処理、磁気データ蓄積処理、磁気データ分布判定処理、中心点算出処理、歪判定処理、及び中心点出力処理を行うが、本発明はこれに限定されるものでは無く、これらのうちの一部について行うものであってもよい。例えば、中心点導出部200は、歪判定処理を行わず、中心点算出処理において算出された中心点cOの座標を、中心点出力処理において出力してもよい。また、中心点導出部200は、磁気データ分布判定処理を行わずに、中心点算出処理を実行してもよい。この場合、中心点導出部200は、磁気センサのオフセットqOFFが急激且つ大きな変化をした場合であっても、変化後の磁気センサのオフセットqOFFの値を容易に算出することができる。3変形例3上述した実施形態では、カルマンフィルタ演算部120は、遅延部121において、状態ベクトルxのうち磁気センサのオフセット推定値qOを表す要素を、中心点導出部200から出力された中心点cOの座標によって上書きするものであるが、本発明はこのような形態に限定されるものでは無く、遅延部121以外において、状態ベクトルxの中心点cOの座標による上書きを行ってもよい。例えば、シグマポイント生成部122において、2a1個のシグマポイントχkを生成する際に、式14及び15の右辺に現れる更新後の状態ベクトルxkのうち磁気センサのオフセット推定値qOを表す要素を、中心点cOの座標によって上書きしてもよい。また、状態遷移モデル部123の演算で用いる状態遷移モデルにおいて、式27の右辺のうち磁気センサのオフセット推定値qOK-1を表す要素を、中心点cOの座標によって上書きしてもよい。さらには、加算器129において算出される、更新後の状態ベクトルxkのうち磁気センサのオフセット推定値qOを表す要素を、中心点cOの座標によって上書きしてもよい。4変形例4上述した実施形態では、カルマンフィルタ演算部120は、非線形カルマンフィルタであるシグマポイントカルマンフィルタによる演算を実行しているが、本発明はこれに限定されるのではなく、カルマンフィルタ演算部120は、拡張カルマンフィルタなどの、シグマポイントカルマンフィルタ以外の非線形カルマンフィルタによる演算を実行するものであってもよい。5変形例5上述した実施形態では、状態ベクトルxを構成する要素して、携帯機器1の姿勢μ、地磁気の強さr、地磁気の伏角φ、携帯機器1の角速度ω、角速度センサのオフセット推定値gO、及び磁気センサのオフセット推定値qOを採用し、これら6個の状態変数をカルマンフィルタの演算における推定の対象としたが、本発明はこれに限定されるものでは無い。例えば、状態ベクトルxを、これら6個の状態変数のうちの一部により構成し、その一部について推定するものであってもよい。6変形例6上述した実施形態では、観測値ベクトルykを、磁気センサの出力である観測値q、加速度センサの出力である観測値a、及び角速度センサの出力である観測値gを要素とするベクトルとしたが、本発明はこれに限定されるものでは無く、このうちの一部のみを利用して観測値ベクトルを生成してもよい。"