"磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出する3次元磁気センサと、前記3次元磁気センサから、3軸の座標系において表現される3次元のベクトルデータとして順次出力される磁気データを蓄積する蓄積手段と、前記蓄積手段に蓄積された前記複数の磁気データで示される複数の座標の分布の3次元的な広がりの程度を示す分散評価値を算出する、分布指標算出手段と、前記複数の磁気データの各々が示す3軸の座標が、歪判定用球面を歪ませた形状の立体の表面近傍に確率的に分布すると仮定して、前記複数の磁気データと、前記分散評価値とに基づいて、前記立体と、前記歪判定用球面との形状の相違の程度を示す歪評価値を算出し、前記歪評価値が歪許容値以下であるか否かを判定する、歪判定手段と、前記歪判定手段の判定結果が肯定である場合に、前記磁気データのうち前記部品の発生する磁界の成分を表す3軸の座標であるオフセットを更新する、オフセット更新手段と、を備えることを特徴とする、地磁気測定装置。前記地磁気測定装置は、前記複数の磁気データの各々が示す3軸の座標が、中心点算出用球面の近傍に確率的に分布すると仮定して、前記複数の磁気データに基づいて前記中心点算出用球面の中心点を算出する、中心点算出手段をさらに備え、前記オフセット更新手段は、前記歪判定手段の判定結果が肯定である場合に、前記オフセットを、旧オフセットから前記中心点算出用球面の中心点に更新する、ことを特徴とする、請求項1に記載の地磁気測定装置。前記歪判定手段の判定結果が否定である場合に、利用者に前記機器の位置を変化させずに前記機器の姿勢を変化させるように促す手段を備えることを特徴とする、請求項1または2に記載の地磁気測定装置。前記分散評価値は、前記複数の磁気データの分散を表す共分散行列の最小固有値であることを特徴とする、請求項1乃至3のうちいずれか1項に記載の地磁気測定装置。前記磁気データを前記中心点算出用球面の中心点を原点とする座標系で表した3軸のベクトルと、当該3軸のベクトルを対称行列である歪評価行列により変換したベクトルとの内積を、前記複数の磁気データの各々について算出し、複数の算出結果を要素とするベクトルを歪誤差ベクトルとし、前記複数の磁気データの各々で示される3軸の座標により特定されるそれぞれの位置と前記歪判定用球面との誤差を表すベクトルを、歪判定用球面誤差ベクトルとし、前記歪誤差ベクトルと、前記歪判定用球面誤差ベクトルとの和を立体誤差ベクトルとし、前記共分散行列の前記最小固有値に対応する正規化された固有ベクトルと、前記旧オフセットから前記歪判定用球面の中心点を示す歪判定用更新ベクトルとの内積に、第1係数を乗じた値を、第1誤差とし、前記歪評価行列と、前記固有ベクトルとの積に、第2係数を乗じた値を第2誤差とし、前記立体誤差ベクトル、前記第1誤差、及び前記第2誤差を要素とするベクトルを、歪評価ベクトルとしたとき、前記歪評価値は、前記歪評価行列の各々の成分と、前記歪判定用球面の中心点を表す3軸のベクトルの各々の要素とを変数とし、前記歪評価ベクトルの大きさを表す歪評価関数を最小化するときの、前記歪評価行列のノルムであることを特徴とする、請求項4に記載の地磁気測定装置。前記複数の磁気データの各々で示される3軸の座標として特定されるそれぞれの位置が、前記中心点算出用球面近傍に確率的に分布すると仮定した場合の、前記複数の磁気データで示される3軸の座標により特定されるそれぞれの位置と前記中心点算出用球面との誤差を表すベクトルを、中心点算出用球面誤差ベクトルとし、前記共分散行列の前記最小固有値に対応する正規化された固有ベクトルと、前記旧オフセットから前記中心点算出用球面の中心点を示す中心点算出用更新ベクトルとの内積に、第1係数を乗じた値を、第3誤差とし、前記中心点算出用球面の中心点は、中心点算出関数の値を最小化する3次元のベクトルで与えられ、前記中心点算出関数は、前記3次元のベクトルを変数とし、且つ、前記中心点算出用球面誤差ベクトルと前記第3誤差とを要素とする中心点算出ベクトルの大きさを表すことを特徴とする、請求項5に記載の地磁気測定装置。前記歪評価関数は、下記のgSD2eで表されることを特徴とする請求項5または6に記載の地磁気測定装置。前記中心点算出関数は、下記のgSxで表されることを特徴とする請求項6または7に記載の地磁気測定装置。前記第1係数は、前記最小固有値が閾値よりも小さい値である場合には、正の第1の値を有し、前記最小固有値が前記閾値以上の場合には、0または前記第1の値よりも小さい正の第2の値を有し、前記第2係数は、前記最小固有値が前記閾値よりも小さい値である場合には、正の第3の値を有し、前記最小固有値が前記閾値以上の場合には、0または前記第3の値よりも小さい正の第4の値を有する、ことを特徴とする、請求項5乃至8のうちいずれか1項に記載の地磁気測定装置。kα及びkβを正の定数としたとき、前記第1係数は、以下のαで表され、前記第2係数は、以下のβで表されることを特徴とする、請求項5乃至8のうちいずれか1項に記載の地磁気測定装置。前記磁気データを前記中心点算出用球面の中心点を原点とする座標系で表した3軸のベクトルと、当該3軸のベクトルを対称行列である歪評価行列により変換したベクトルとの内積を、前記複数の磁気データの各々について算出し、複数の算出結果を要素とするベクトルを歪誤差ベクトルとし、前記複数の磁気データの各々で示される3軸の座標として特定されるそれぞれの位置と前記歪判定用球面との誤差を表すベクトルを、歪判定用球面誤差ベクトルとし、前記歪誤差ベクトルと、前記歪判定用球面誤差ベクトルとの和を、立体誤差ベクトルとし、前記歪評価行列の各々の成分と、前記歪判定用球面の中心点を表す3軸のベクトルの各々の要素とを変数とし、前記立体誤差ベクトルの大きさを表す関数を、立体歪評価関数とし、前記共分散行列の前記最小固有値に対応する正規化された固有ベクトルと、前記旧オフセットから前記歪判定用球面の中心点を示す歪判定用更新ベクトルとが直交する条件を、第1の制約条件とし、前記歪評価行列と、前記固有ベクトルとの積が、零ベクトルとなる条件を、第2の制約条件としたとき、前記最小固有値が閾値よりも小さい値である場合、前記歪評価値は、前記第1の制約条件及び前記第2の制約条件を充足したうえで前記立体歪評価関数を最小化するときの、前記歪評価行列のノルムであることを特徴とする、請求項4に記載の地磁気測定装置。前記立体歪評価関数は、以下のgSDeで表され、前記第1の制約条件は、以下の式aで表され、前記第2の制約条件は、以下の式bで表されることを特徴とする、請求項11に記載の地磁気測定装置。前記歪評価行列Eは、以下の式cの解として求められるベクトルe=[e11、e22、e12、e33、e23、e13、eXT]Tのうち、第1行目から第6行目までの要素を成分とすることを特徴とする、請求項12に記載の地磁気測定装置。磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出する3次元磁気センサと、前記3次元磁気センサから、3軸の座標系において表現される3次元のベクトルデータとして順次出力される磁気データを蓄積する蓄積手段と、を備える地磁気測定装置に用いられるオフセット更新方法であって、前記蓄積手段に蓄積された前記複数の磁気データで示される複数の座標の分布の、3次元的な広がりの程度を示す分散評価値を算出し、前記複数の磁気データの各々が示す3軸の座標が、歪判定用球面を歪ませた形状の立体の表面近傍に確率的に分布すると仮定して、前記複数の磁気データと、前記分散評価値とに基づいて、前記立体と、前記歪判定用球面との形状の相違の程度を示す歪評価値を算出し、前記歪評価値が歪許容値以下であるか否かを判定し、前記判定の結果が肯定である場合に、前記磁気データのうち前記部品の発生する磁界の成分を表す3軸の座標であるオフセットを更新する、ことを特徴とする、オフセット更新方法。磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出する3次元磁気センサと、前記3次元磁気センサから、3軸の座標系において表現される3次元のベクトルデータとして順次出力される磁気データを蓄積する蓄積手段と、を備える地磁気測定装置に用いられるオフセット更新プログラムであって、前記蓄積手段に蓄積された前記複数の磁気データで示される複数の座標の分布の、3次元的な広がりの程度を示す分散評価値を算出する処理と、前記複数の磁気データの各々が示す3軸の座標が、歪判定用球面を歪ませた形状の立体の表面近傍に確率的に分布すると仮定して、前記複数の磁気データと、前記分散評価値とに基づいて、前記立体と、前記歪判定用球面との形状の相違の程度を示す歪評価値を算出し、前記歪評価値が歪許容値以下であるか否かを判定する処理と、前記判定の結果が肯定である場合に、前記磁気データのうち前記部品の発生する磁界の成分を表す3軸の座標であるオフセットを更新する処理とを、コンピュータに実行させるオフセット更新プログラム。ところで、3次元磁気センサが搭載される携帯電話等の機器の近くに、スピーカやパソコン等の磁界を発生させる物体が存在する場合、3次元磁気センサは、前述した地磁気及び内部磁界の他に、機器の外部の物体が発する外部磁界も検出する。一般的に、外部磁界は、外部磁界を発生させる物体と3次元磁気センサとの相対的位置関係等によって、向き及び大きさが変化する不均一な磁界である。従来の方法では、このような不均一な外部磁界の存在の有無や、不均一な外部磁界による影響の大きさを知ることは出来なかった。従って、従来の方法では、不均一な外部磁界が存在する環境下で算出された内部磁界のベクトルが大きな誤差を持つ場合であっても、外部磁界が存在しない環境下で計算された正しい値であると判断し、これをオフセットとして採用してしまうため、オフセットに大きな誤差を含むことがあるという問題が存在した。また、特許文献1に開示されているような、機器の姿勢変化が不十分な場合に、3次元磁気センサを姿勢変化させつつ取得した複数の磁気データと、過去に算出したオフセットとを用いることによりオフセットを算出する方法も、不均一な外部磁界の存在を考慮するものではなく、不均一な外部磁界が存在する環境では、外部磁界の影響を受けた誤差を含むオフセットを、正確なオフセットであると判断してしまうという問題が存在した。そこで、本発明は上述した事情を鑑み、機器の姿勢変化が不十分な場合にも、外部磁界の影響の程度を評価することより、算出したオフセット値が誤差を含まない正確な値であるかを判断し、算出したオフセットの値の採否を決定することを解決課題とする。以下、本発明について説明する。なお、本発明の理解を容易にするために本実施形態、変形例、及び添付図面の参照符号を括弧書きにて付記するが、それにより本発明が本実施形態に限定されるものではない。上述した課題を解決するため、本発明に係る地磁気測定装置は、磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出する3次元磁気センサと、前記3次元磁気センサから、3軸の座標系において表現される3次元のベクトルデータとして順次出力される磁気データqiを蓄積する蓄積手段と、前記蓄積手段に蓄積された前記複数の磁気データq1~qNで示される複数の座標の分布の3次元的な広がりの程度を示す分散評価値を算出する、分布指標算出手段と、前記複数の磁気データの各々が示す3軸の座標が、歪判定用球面S2を歪ませた形状の立体SDの表面近傍に確率的に分布すると仮定して、前記複数の磁気データと、前記分散評価値とに基づいて、前記立体SDと、前記歪判定用球面S2との形状の相違の程度を示す歪評価値gDとの内積を、前記複数の磁気データの各々について算出し、複数の算出結果を要素とするベクトルを歪誤差ベクトルkとし、前記複数の磁気データの各々で示される3軸の座標により特定されるそれぞれの位置と前記歪判定用球面S2との誤差を表すベクトルを、歪判定用球面誤差ベクトルδS2とし、前記歪誤差ベクトルkと、前記歪判定用球面誤差ベクトルδS2との和を立体誤差ベクトルδSDとし、前記共分散行列の前記最小固有値λ3に対応する正規化された固有ベクトルu3と、前記旧オフセットから前記歪判定用球面の中心点x02を示す歪判定用更新ベクトルk2との内積に、第1係数αを乗じた値を、第1誤差δU2とし、前記歪評価行列と、前記固有ベクトルu3との積に、第2係数βを乗じた値を第2誤差δEとし、前記立体誤差ベクトル、前記第1誤差、及び前記第2誤差を要素とするベクトルを、歪評価ベクトルδSD2としたとき、前記歪評価値gDEは、前記歪評価行列の各々の成分e11~e33と、前記歪判定用球面の中心点x02を表す3軸のベクトルの各々の要素とを変数とし、前記歪評価ベクトルの大きさを表す歪評価関数gSD2を最小化するときの、前記歪評価行列のノルムであることを特徴とする、ことが好ましい。この発明によれば、立体誤差ベクトルを要素として含む歪評価ベクトルの大きさを最小化することで、歪評価行列、及び歪判定用球面の中心点を算出する。立体誤差ベクトルは、歪判定用球面誤差ベクトルと、歪誤差ベクトルとの和で与えられる。ここで、歪判定用球面誤差ベクトルは、複数の磁気データで示される座標と歪判定用球面との誤差を、任意の値を取り得る制約の無い変数で表した値を各要素とするベクトルである。従って、歪判定用球面誤差ベクトルの各要素、すなわち、複数の磁気データに対応する誤差の各々は、互いに独立なものである。その結果、複数の磁気データで示される座標と歪判定用球面との誤差を最小化するように歪判定用球面を定めた場合、歪判定用球面誤差ベクトルにより表される複数の誤差は、全体として対称性を有し、方向依存性の無い、ホワイトノイズとなる。一方、歪誤差ベクトルは、中心点算出用球面の中心点から見た磁気データの示す座標を表すベクトルと、当該ベクトルを歪評価行列により変換したベクトルとの内積を、複数の磁気データの各々について並べたベクトルである。歪評価行列は対称行列であるため、歪誤差ベクトルの各要素は、対称行列である歪評価行列を係数行列とし、中心点から見た磁気データの座標を表すベクトルの3つの要素の各々を変数とする、3変数二次形式となる。すなわち、歪誤差ベクトルは、複数の磁気データで示される座標と歪判定用球面との誤差の各々を、二次形式で表される同一の関数に基づく曲面上に存在するという制約の下で表現するベクトルである。このような歪誤差ベクトルを用いて、複数の磁気データと歪判定用球面との誤差を表現することにより、ホワイトノイズとは別に、二次形式の関数に基づく曲面で表される誤差すなわち、球面からの歪を表現することが可能となり、これまで一律にホワイトノイズとみなしていた誤差から、ホワイトノイズとは異なる球面からの歪による誤差を分離することが可能となる。歪判定用球面誤差ベクトルと歪誤差ベクトルとを加算して得られる立体誤差ベクトルは、複数の磁気データを表面近傍に有し、球面とは異なる形状の立体を表現する。そして、立体誤差ベクトルのうち、歪誤差ベクトルの大きさ評価することで、立体の形状が、歪判定用球面の形状からどの程度の差異歪を有するかつまり、ホワイトノイズとは異なる球面からの歪による誤差の大きさを、評価することが可能となる。その結果、立体と歪判定用球面とが実質的に同一の形状を有するほど近似していると評価される場合には、複数の磁気データの示す座標が歪判定用球面近傍にも分布すると看做すことができる。前述の通り、複数の磁気データの示す座標が、球面の形状となるように分布する場合、複数の磁気データの示す座標は、中心点算出用球面近傍にも分布すると看做すことができる。従って、立体と歪判定用球面とが実質的に同一の形状を有するほど近似していると評価される場合には、中心点算出用球面の中心点をオフセットとして採用する。一方、立体が歪判定用球面とは大きく異なる歪んだ形状を有すると評価される場合には、複数の磁気データで示される座標が歪判定用球面近傍に分布すると看做すことはできないため、中心点算出用球面の中心点をオフセットとして採用することを防止する。すなわち、上述した地磁気測定装置は、中心点算出手段によって算出された中心点のうち、不均一な外部磁界の影響を受けていない複数の磁気データより算出された、適切な中心点のみをオフセットとして採用することが可能であるという利点を有する。この結果、正確な地磁気を検出することが可能となる。本発明は、立体と球面との形状の相違の程度を示す歪評価値を、歪評価行列のノルムとして定めている。歪評価行列は、3行3列の対称行列であり、3つの固有値と、互いに直交する3つの固有ベクトルとを有する。歪評価行列のノルムは、歪評価行列が有する3つの固有値のうち、絶対値が最大となる固有値に基づいて定められる。ところで、歪評価関数を最小化することで算出される歪評価行列は、複数の磁気データの示す座標を表面近傍に有する立体と歪判定用球面との形状の相違を表す歪誤差ベクトルを定める行列である。ここで、歪誤差ベクトルは、複数の磁気データの示す座標と歪判定用球面との誤差のうちホワイトノイズとして吸収しきれない誤差、すなわち、立体と歪判定用球面との形状の相違を表すベクトルである。つまり、歪誤差ベクトルを構成する複数の要素のうち、ある要素の絶対値は、当該要素に対応する磁気データの座標と歪判定用球面との誤差が大きい場合には大きな値に定められ、当該要素に対応する磁気データの座標と歪判定用球面との誤差が小さい場合には小さな値に定められる。従って、歪評価行列の各成分は、歪誤差ベクトルを構成する複数の要素が全体として、複数の磁気データの座標と歪判定用球面との誤差の大きさを正確に表現することができるように定められる。ところで、歪誤差ベクトルの各要素は、中心点算出用球面の中心点から見た磁気データの示す座標を表すベクトルと、当該ベクトルを歪評価行列により変換したベクトルとの内積として与えられる。つまり、歪誤差ベクトルを構成する複数の要素のうち、ある要素の絶対値は、当該要素に対応する磁気データの示す座標を中心点算出用球面の中心点から表したベクトルと、歪評価行列が有する3つの固有値のうち絶対値の大きさが最大となる固有値に対応する固有ベクトルとが平行となる場合には、大きな値となる。従って、歪評価行列が有する3つの固有値のうち絶対値が最大となる固有値に対応する固有ベクトルが、歪判定用球面からの距離が大きな位置に存在する磁気データが多く存在する領域を示すように、歪評価行列の各成分を定める場合、歪誤差ベクトルは、複数の磁気データの示す座標と歪判定用球面との誤差の大きさを正確に表現できることになる。このように、歪評価行列の各成分は、絶対値の大きさが最大となる固有値に対応する固有ベクトルの方向と、歪判定用球面からの誤差が大きな磁気データが多く存在する領域を中心点算出用球面の中心点を起点として表したベクトルの方向とが、近づくように定められる。そして、歪評価行列が有する3つの固有値のうち、絶対値の大きさが最大となる固有値は、立体の表面近傍に存在する磁気データの座標と球面との誤差が大きい場合に、その誤差の大きさを表現する値となる。本発明は、立体と球面との形状の相違の程度を示す歪評価値を、歪評価行列の有する3つの固有値のうち絶対値が最大となる固有値の値に基づいて定めている。従って、歪評価値は、立体と球面との誤差の程度の大きさを評価することが可能となる。前述の通り、機器の動かし方が不十分で複数の磁気データの示す座標がある平面上に2次元的に分布する場合には、立体の3次元的な形状を把握することは出来ない。この場合には、立体を平面で切断した切断面を表す曲線と、球面を平面で切断した切断面を表す円周との形状の相違の程度を判断することで、外部磁界の影響の大きさを判断することが必要となる。この発明によれば、第2誤差を要素として含む歪評価ベクトルの大きさを最小化することで、歪評価行列を算出する。歪評価ベクトルの大きさを最小化して0に近づける場合、第2誤差も最小化されて0に近づく。第2誤差は、共分散行列が有する3つの固有値のうち最小の固有値最小固有値に対応する固有ベクトルを歪評価行列により変換したベクトルに対して、第2係数を乗じたベクトルである。従って、第2誤差を最小化して0に近づける場合、共分散行列の最小固有値に対応する固有ベクトルを歪評価行列により変換したベクトルの大きさ、または、第2係数のうち、少なくともいずれか一方が最小化されて0に近づく。ここで、共分散行列の最小固有値に対応する固有ベクトルを歪評価行列により変換したベクトルの大きさが最小化されて0に近づく場合を検討する。共分散行列の最小固有値に対応する固有ベクトルを歪評価行列により変換したベクトルの大きさが最小化されて0に近づくという条件の下で、歪評価行列を算出する場合、歪評価行列の有する3つの固有値のうち、絶対値が最小の固有値は0または0に近い値となる。同時に、歪評価行列の有する3つの固有値のうち、絶対値が最小の固有値に対応する固有ベクトルは、共分散行列の最小固有値に対応する固有ベクトルと平行になる。複数の磁気データがある平面上に2次元的に分布する場合、共分散行列の最小固有値に対応する固有ベクトルは、磁気データが分布する平面に直交する。従って、歪評価行列の有する3つの固有ベクトルのうち、絶対値が最小の固有値に対応する固有ベクトルも、磁気データが分布する平面に直交する。つまり、共分散行列の最小固有値に対応する固有ベクトルを歪評価行列により変換したベクトルの大きさが最小化されて0に近づくという条件の下で算出された歪評価行列は、磁気データが分布する平面に直交する固有ベクトル方向については、立体と歪判定用球面との形状の相違は無いものと看做し、立体と歪判定用球面との形状の相違を評価しない行列である。一方、歪評価行列の有する3つの固有ベクトルは互いに直交するため、絶対値が最小の固有値に対応する固有ベクトル以外の2つの固有ベクトルは、磁気データが分布する平面と平行になる。そして、磁気データが分布する平面に平行な2つの固有ベクトルに対応する固有値は、0以外の値を取ることができる。従って、共分散行列の最小固有値に対応する固有ベクトルを歪評価行列により変換したベクトルの大きさが最小化されて0に近づくという条件の下で算出された歪評価行列は、磁気データが分布する平面に平行な2つの固有ベクトル方向について、立体と歪判定用球面との形状の相違の程度を評価することができる行列である。このように、複数の磁気データがある平面上に2次元的に分布する場合には、共分散行列の最小固有値に対応する固有ベクトルを歪評価行列により変換したベクトルの大きさを最小化するという条件の下で算出された歪評価行列を用いることで、磁気データが分布する平面と平行な方向については、立体と球面との形状の相違の程度を評価し、且つ、磁気データが分布する平面と垂直な方向については、立体と球面との形状の相違の程度を評価しないことが可能となる。なお、磁気データの示す座標が3次元的に分布する場合には、共分散行列の最小固有値に対応する固有ベクトル方向についても、立体と歪判定用球面との形状の相違の程度を評価することが必要となる。つまり、この場合、歪評価行列は、共分散行列の最小固有値に対応する固有ベクトル方向についても、立体と歪判定用球面との形状の相違を評価できる行列である必要がある。この発明によれば、第2誤差は、共分散行列の最小固有値に対応する固有ベクトルを歪評価行列により変換したベクトルに第2係数を乗じたベクトルである。よって、第2係数を0または0に近い小さな値に設定することで、共分散行列の最小固有値に対応する固有ベクトルを歪評価行列により変換したベクトルの大きさを最小化することなく、第2誤差を最小化することができる。この場合、歪評価行列は、絶対値が最小の固有値も0以外の値を取ることができる。従って、歪評価行列は、3次元空間内の全ての方向について、立体と歪判定用球面との形状の相違の程度を表現することができる。このように、この発明は、第2係数を有する第2誤差を最小化することで歪評価行列を算出する。そして、複数の磁気データの示す座標の分布の3次元的な広がりの程度に応じて、第2係数の値を設定することで、共分散行列の最小固有値に対応する固有ベクトルを歪評価行列により変換したベクトルの大きさを最小化するという条件を、緩和することが可能となる。すなわち、本発明の歪評価行列を用いることで、複数の磁気データの示す座標が3次元的に分布する場合には、立体と球面との形状の相違の程度を3次元的に評価し、複数の磁気データの示す座標がある平面上に2次元的に分布する場合には、立体と球面との形状の相違の程度を、2次元的に当該ある平面内で評価することが可能となる。これにより、磁気データの示す座標の分布が、3次元的であるか2次元的であるかに関わらず、不均一な外部磁界の影響の大きさを評価することが可能となり、不正確なオフセットを採用することを防止できる。また、上述した地磁気測定装置において、前記複数の磁気データの各々で示される3軸の座標として特定されるそれぞれの位置が、前記中心点算出用球面S近傍に確率的に分布すると仮定した場合の、前記複数の磁気データで示される3軸の座標により特定されるそれぞれの位置と前記中心点算出用球面Sとの誤差を表すベクトルを、中心点算出用球面誤差ベクトルδSとし、前記共分散行列の前記最小固有値λ3に対応する正規化された固有ベクトルu3と、前記旧オフセットから前記中心点算出用球面の中心点x0を示す中心点算出用更新ベクトルkとの内積に、第1係数αを乗じた値を、第3誤差δUとし、前記中心点算出用球面の中心点x0は、中心点算出関数gSの値を最小化する3次元のベクトルxで与えられ、前記中心点算出関数は、前記3次元のベクトルを変数とし、且つ、前記中心点算出用球面誤差ベクトルδSと前記第3誤差δUとを要素とする中心点算出ベクトルδSUの大きさを表すことを特徴とする、ことが好ましい。複数の磁気データの示す座標がある平面上に2次元的に分布する場合に、複数の磁気データの示す座標が3次元的に分布すると仮定して算出された中心点の座標は、平面に垂直な方向について不正確である。従って、複数の磁気データの示す座標が2次元的に分布する場合には、複数の磁気データの示す座標が3次元的に分布すると仮定して算出された球面の中心点の座標が有する情報のうち、不正確な情報である平面に垂直な方向共分散行列の最小固有値に対応する固有ベクトル方向についての情報を排除して、中心点を算出する必要がある。この発明によれば、第3誤差を要素として含む歪評価ベクトルの大きさを最小化することで、中心点算出用球面の中心点を算出する。第3誤差は、共分散行列の最小固有値に対応する固有ベクトルと、旧オフセットから中心点算出用球面の中心点の座標を表す中心点算出用更新ベクトルとの内積に、第1係数を乗じた値である。従って、第3誤差を最小化して0に近づける場合、共分散行列の最小固有値に対応する固有ベクトルと中心点算出用更新ベクトルとの内積、または、第1係数のうち、少なくともいずれか一方が最小化されて0に近づく。共分散行列の最小固有値に対応する固有ベクトルと中心点算出用更新ベクトルとの内積が最小化され0に近づく場合、共分散行列の最小固有値に対応する固有ベクトルと中心点算出用更新ベクトルとは直交する。つまり、中心点算出用更新ベクトルは、共分散行列の最小固有値に対応する固有ベクトル方向の成分を含まないベクトルになる。中心点算出用球面の中心点は、旧オフセットと中心点算出用更新ベクトルとの和であるため、不正確な情報共分散行列の最小固有値に対応する固有ベクトル方向の情報を利用せずに算出される。その結果、複数の磁気データが2次元的に分布する場合に、正確な情報のみに基づいて中心点を算出することが可能となる。一方、磁気データの示す座標が3次元的に分布する場合には、第1係数を0または0に近い小さな値に設定することで、共分散行列の最小固有値に対応する固有ベクトルと中心点算出用更新ベクトルとの内積を最小化することなく、第3誤差を最小化することができる。すなわち、第1係数によって、共分散行列の最小固有値に対応する固有ベクトルと中心点算出用更新ベクトルとが直交するという条件を緩和することができる。これにより、複数の磁気データの示す座標が3次元的な広がりを有する場合には、共分散行列の最小固有値に対応する固有ベクトル方向の情報を利用して中心点算出用球面の中心点の座標を算出することが可能となる。また、上述した地磁気測定装置において、前記歪評価関数は、下記のgSD2eで表されることが好ましい。この発明によれば、歪評価関数を最小化するときの歪評価行列のノルムである歪評価値を簡単な計算により算出することができる。従って、上述した地磁気測定装置は、処理速度の向上という利点を有する。また、上述した地磁気測定装置において、前記中心点算出関数は、下記のgSxで表されることが好ましい。この発明によれば、簡単な計算により中心点を算出することができる。従って、上述した地磁気測定装置は、処理速度の向上という利点を有する。また、上述した地磁気測定装置において、前記第1係数αは、前記最小固有値λ3が閾値λ0よりも小さい値である場合には、正の第1の値を有し、前記最小固有値が前記閾値以上の場合には、0または前記第1の値よりも小さい正の第2の値を有し、前記第2係数βは、前記最小固有値が前記閾値よりも小さい値である場合には、正の第3の値を有し、前記最小固有値が前記閾値以上の場合には、0または前記第3の値よりも小さい正の第4の値を有する、ことを特徴とする、ことが好ましい。この発明によれば、共分散行列の最小固有値が閾値よりも小さな値を有し、複数の磁気データの示す座標がある平面上に2次元的に分布すると看做せる場合には、第1係数及び第2係数を正の値とする。第1係数を正の第1の値とすることで、共分散行列の最小固有値に対応する固有ベクトルと中心点算出用更新ベクトルとが直交するという条件の下、正確な情報のみに基づいて、中心点算出用球面の中心点を算出することが可能となる。また、第1係数を正の第1の値とすることで、共分散行列の最小固有値に対応する固有ベクトルと歪判定用更新ベクトルとが直交するという条件の下、正確な情報のみに基づいて、歪判定用球面の中心点を算出することが可能となる。さらに、第2係数を正の第3の値とした場合、共分散行列の最小固有値に対応する固有ベクトルを歪評価行列により変換したベクトルの大きさが0に近い値となるという条件の下で歪評価行列が算出される。このような歪評価行列によって、複数の磁気データの示す座標が分布する平面に平行な方向についてのみ立体と球面の形状の相違の程度を評価することができる。一方、共分散行列の最小固有値が閾値以上となり、複数の磁気データの示す座標が3次元的な広がりを有して分布する場合には、第1係数及び第2係数を0または0に近い小さな値とする。第1係数を0または0に近い値とすることで、旧オフセットの情報を用いずに、複数の磁気データのみに基づいて中心点算出用球面の中心点及び歪判定用球面の中心点を算出することが可能となる。また、第2係数を0または0に近い値として算出された歪評価行列により、立体と球面との3次元的な形状の相違の程度を、評価することが可能になる。これにより、複数の磁気データの示す座標の分布が3次元的な場合及び2次元的な場合の双方において、外部磁界の影響を受けない適切な中心点をオフセットに採用することが可能となる。また、上述した地磁気測定装置において、kα及びkβを正の定数としたとき、前記第1係数αは、以下のαで表され、前記第2係数βは、以下のβで表されることが好ましい。この発明によれば、磁気データの示す座標の分布の3次元的な広がりの程度に応じて、歪評価行列と共分散行列の最小固有値に対応する固有ベクトルとの積が0に近い値となるという条件を緩和するため、複数の磁気データの示す座標の分布の広がりの程度に応じて適切に立体と球面との形状の相違の程度を評価することが可能となる。また、複数の磁気データの示す座標の分布の3次元的な広がりの程度に応じて、中心点算出用更新ベクトルと共分散行列の最小固有値に対応する固有ベクトルとが直交するという条件、及び、歪判定用更新ベクトルと共分散行列の最小固有値に対応する固有ベクトルとが直交するという条件を緩和するため、複数の磁気データの示す座標の3次元的な分布の広がりの程度に応じて、正確な中心点算出用球面の中心点及び歪判定用球面の中心点を算出することが可能となる。また、上述した地磁気測定装置において、前記磁気データを前記中心点算出用球面の中心点x0を原点とする座標系で表した3軸のベクトルと、当該3軸のベクトルを対称行列である歪評価行列Eにより変換したベクトルとの内積を、前記複数の磁気データの各々について算出し、複数の算出結果を要素とするベクトルを歪誤差ベクトルkとし、前記複数の磁気データの各々で示される3軸の座標として特定されるそれぞれの位置と前記歪判定用球面S2との誤差を表すベクトルを、歪判定用球面誤差ベクトルδS2とし、前記歪誤差ベクトルkと、前記歪判定用球面誤差ベクトルδS2との和を、立体誤差ベクトルδSDとし、前記歪評価行列の各々の成分と、前記歪判定用球面の中心点x02を表す3軸のベクトルの各々の要素とを変数とし、前記立体誤差ベクトルの大きさを表す関数を、立体歪評価関数gSDとし、前記共分散行列の前記最小固有値λ3に対応する正規化された固有ベクトルu3と、前記旧オフセットから前記歪判定用球面の中心点x02を示す歪判定用更新ベクトルk2とが直交する条件を、第1の制約条件とし、前記歪評価行列と、前記固有ベクトルとの積が、零ベクトルとなる条件を、第2の制約条件としたとき、前記最小固有値λ3が閾値λ0よりも小さい値である場合、前記歪評価値gDEは、前記第1の制約条件及び前記第2の制約条件を充足したうえで前記立体歪評価関数gSDを最小化するときの、前記歪評価行列のノルムであることを特徴とする、ことが好ましい。この発明によれば、最小固有値が閾値よりも小さな値を有し、複数の磁気データの示す座標がある平面上に2次元的に分布すると看做すことができる場合には、共分散行列の最小固有値に対応する固有ベクトルと歪判定用更新ベクトルとが直交するという第1の制約条件と、共分散行列の最小固有値に対応する固有ベクトルを歪評価行列により変換したベクトルの大きさが0になるという第2の制約条件とを厳格に満たしつつ、歪評価行列及び歪判定用球面の中心点を算出し、立体と歪判定用球面との形状の相違の程度を評価する。従って、共分散行列の最小固有値に対応する固有ベクトル方向の不正確な情報を完全に排除して算出される歪判定用球面の中心点の情報を用いて、平面と完全に平行な方向についてのみ立体と球面の形状の相違の程度を評価することが可能となる。これにより、複数の磁気データの示す座標が2次元的に分布する場合に、立体と球面の形状の相違の程度をより正確に評価することが可能となる。また、上述した地磁気測定装置において、前記立体歪評価関数は、以下のgSDeで表され、前記第1の制約条件は、以下の式aで表され、前記第2の制約条件は、以下の式bで表されることを特徴とことが好ましい。この発明によれば、歪評価関数を最小化するときの歪評価行列のノルムである歪評価値を、立体誤差ベクトル、第1の制約条件、及び第2の制約条件に基づいて簡単な計算により算出することができる。従って、上述した地磁気測定装置は、処理速度の向上という利点を有する。また、上述した地磁気測定装置において、前記歪評価行列Eは、以下の式cの解として求められるベクトルe=[e11、e22、e12、e33、e23、e13、eXT]Tのうち、第1行目から第6行目までの要素を成分とすることが好ましい。この発明によれば、簡単な計算により歪評価行列を算出することができる。従って、上述した地磁気測定装置は、処理速度の向上という利点を有する。次に、本発明に係る地磁気測定装置に用いられるオフセット更新方法は、磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出する3次元磁気センサと、前記3次元磁気センサから、3軸の座標系において表現される3次元のベクトルデータとして順次出力される磁気データを蓄積する蓄積手段と、を備える地磁気測定装置に用いられるオフセット更新方法であって、前記蓄積手段に蓄積された前記複数の磁気データで示される複数の座標の分布の、3次元的な広がりの程度を示す分散評価値を算出し、前記複数の磁気データの各々が示す3軸の座標が、歪判定用球面を歪ませた形状の立体の表面近傍に確率的に分布すると仮定して、前記複数の磁気データと、前記分散評価値とに基づいて、前記立体と、前記歪判定用球面との形状の相違の程度を示す歪評価値を算出し、前記歪評価値が歪許容値以下であるか否かを判定し、前記判定の結果が肯定である場合に、前記磁気データのうち前記部品の発生する磁界の成分を表す3軸の座標であるオフセットを更新する、ことを特徴とする。この発明によれば、複数の磁気データの示す座標が3次元的に分布する場合、及び2次元的に分布する場合の双方について、立体と球面との形状の相違の程度を評価して、外部磁界の影響の受けていない正確なオフセットを採用することができる。次に、本発明に係る地磁気測定装置に用いられるオフセット更新プログラムは、磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出する3次元磁気センサと、前記3次元磁気センサから、3軸の座標系において表現される3次元のベクトルデータとして順次出力される磁気データを蓄積する蓄積手段と、を備える地磁気測定装置に用いられるオフセット更新プログラムであって、前記蓄積手段に蓄積された前記複数の磁気データで示される複数の座標の分布の、3次元的な広がりの程度を示す分散評価値を算出する処理と、前記複数の磁気データの各々が示す3軸の座標が、歪判定用球面を歪ませた形状の立体の表面近傍に確率的に分布すると仮定して、前記複数の磁気データと、前記分散評価値とに基づいて、前記立体と、前記歪判定用球面との形状の相違の程度を示す歪評価値を算出し、前記歪評価値が歪許容値以下であるか否かを判定する処理と、前記判定の結果が肯定である場合に、前記磁気データのうち前記部品の発生する磁界の成分を表す3軸の座標であるオフセットを更新する処理とを、コンピュータに実行させることを特徴とする。この発明によれば、複数の磁気データの示す座標が3次元的に分布する場合、及び2次元的に分布する場合の双方について、立体と球面との形状の相違の程度を評価して、外部磁界の影響の受けていない正確なオフセットを採用することができる。<A.実施形態>以下、本発明の実施の形態を説明する。[1.3次元磁気センサが検出する磁界の概要]本実施形態では、3次元磁気センサが検出する磁界として、検出対象である地磁気の他に、3次元磁気センサが搭載される機器の部品が発する磁界である内部磁界、及び機器の外部の物体が発する外部磁界の存在を想定している。これらの3種類の磁界は、磁界の向き、大きさ、及び均一性において、各々異なる性質を有するものである。つまり、これらの磁界を3次元磁気センサで検出する場合、3次元磁気センサの位置、姿勢が変化すると、3次元磁気センサから出力される磁界を示すベクトルデータの向き及び大きさも変化するが、この出力データの変化の仕方は、3種類の磁界の各々で差異を有する。以下、図1乃至図4を用いて、本実施形態において想定するこれら3種類の磁界の概要、及び、これらの磁界を3次元磁気センサが検出した場合に出力されるベクトルデータについて説明する。図1は、機器1に搭載された3次元磁気センサ60の周囲に存在する磁界の種類及び性質について説明する図である。3次元磁気センサ60の周囲に存在する磁界には、計測対象である地磁気Bgの他に、3次元磁気センサ60の搭載された機器1の部品が発する内部磁界Bi、及び機器1の周囲に存在する物体2が発する外部磁界Bxが存在する。地磁気Bgは、磁極北に向かう方向を持つベクトルであり、一様な向き及び大きさを有する磁界である。厳密には、地磁気Bgの向き及び大きさは、地域によって異なるが、例えば、異なる都市に移動する等の大きな移動をしない場合には、一様な向き及び大きさを有する。内部磁界Biは、機器1の部品が発する磁界であり、内部磁界Biの向きは、機器1から見て一定の方向を有する。すなわち、機器1上に搭載された3次元磁気センサ60から見た場合、機器1の姿勢を変化させても、内部磁界Biは、一定の方向及び大きさを有する磁界として検知される。例えば、利用者が機器1を手で持って手首を回すと、機器1の姿勢が変化する。このとき、地磁気Bgは磁極北を向いているので、3次元磁気センサ60で検出される地磁気Bgの3軸各々の成分は変化する。一方、3次元磁気センサ60において3次元磁気センサ60と内部磁界を発生する部品との相対的な位置関係は固定であるから、機器1の姿勢を変化させても、内部磁界Biの3軸各々の成分は一定である。外部磁界Bxは、外部の物体2の発する磁界であり、物体2からの距離及び方向によって、磁界の向き及び大きさが変化する不均一な磁界である。すなわち、機器1上に搭載された3次元磁気センサ60から見た場合、外部磁界Bxは、物体2と3次元磁気センサ60との相互位置関係に依存して、方向及び大きさを変化させる磁界として検知される。説明の便宜上、図1に示すような地上座標系ΣG及びセンサ座標系ΣSを導入する。図1に記載された各ベクトルの左上に付された添字Gは、当該ベクトルが地上座標系ΣGにおいて表現されたベクトルであることを意味する。地上座標系ΣGは、地上に固定された座標系であり、地上の任意の一点を原点とし、互いに直交する3つの方向、例えば、東、北、及び鉛直上向きを、それぞれx軸、y軸、及びz軸とする座標系である。ここで、地上座標系ΣGから見たときの3次元磁気センサ60の位置及び姿勢を、それぞれ、位置Ps、及び姿勢θsと表現することにする。センサ座標系ΣSは3次元磁気センサ60に固定された座標系であり、3次元磁気センサ60の有する3つの各々のセンサモジュールの出力する値を、それぞれ、x軸、y軸、z軸上にプロットするように設けられた座標系である。すなわち、3次元磁気センサ60が出力する磁気データは、センサ座標系ΣSのベクトルデータとして表現される。なお、地上座標系ΣGから見た場合、センサ座標系ΣSの原点は位置Psであり、姿勢は姿勢θsである。図2は、地磁気Bg及び内部磁界Biが、3次元磁気センサ60によってどのような値として検知されるかについて示す図である。すなわち、図2は、3次元磁気センサ60の地上座標系ΣGから見た姿勢θsをθ1~θNと変化させつつ磁界を測定した場合に、3次元磁気センサ60が出力するN個の磁気データq1~qNを、センサ座標系ΣSにおいてプロットした図であるNは精度よくオフセットを導出するために必要な磁気データの規定測定回数を表す4以上の自然数。なお、図2では、地磁気Bg及び内部磁界Biのみが存在し、外部磁界Bxは存在しない場合を想定している。また、図2に記載された各ベクトルの左上に付された添字Sは、当該ベクトルがセンサ座標系ΣSにおいて表されたベクトルであることを意味する。内部磁界Biは、センサ座標系ΣSにおいて、3次元磁気センサ60の位置Ps及び姿勢θsに依存しない、一定の方向及び大きさを持つベクトルSBiセンサ座標系ΣSの原点から中心点x0Gへと向かうベクトルとして表現される。一方、地磁気Bgは、センサ座標系ΣSにおいて、大きさは一定であるが、向きは3次元磁気センサ60の姿勢θsに伴い変化する。つまり、地磁気Bgは、センサ座標系ΣSにおいて、3次元磁気センサ60の姿勢θsに依存した向きを有する、一定の大きさのベクトルSBgθsとして表現され、θsを変化させた場合には、中心点x0Gを中心とする球面上に位置することになる。従って、3次元磁気センサ60が順次出力する磁気データq1~qNは、センサ座標系ΣSにおいて、中心点x0Gを中心点とし、地磁気Bgの大きさを半径とする、球面SG上に分布することになる。仮に、内部磁界Biの大きさが0の場合には、中心点x0Gと、センサ座標系ΣSの原点とは一致することになる。なお、3次元磁気センサ60の測定値は、測定誤差を有するため、複数の磁気データq1~qNは、厳密には、球面SGの近傍に確率的に分布する。このような複数の磁気データq1~qNから、球面SGの中心点x0Gの座標を計算することで、内部磁界を表すベクトルSBiを求めることができ、3次元磁気センサ60の出力値である複数の磁気データq1~qNから内部磁界SBiを引き算することにより、地磁気SBgを計算することが可能となる。すなわち、3次元磁気センサ60で検出される磁気が地磁気Bg及び内部磁界Biである場合、球面SGの中心点x0Gを算出することが可能であり、中心点x0Gは3次元磁気センサ60のオフセットcに対応する。図3は、外部磁界Bxが存在する場合に、3次元磁気センサ60の位置PsをP1~PNと変化させると共に、姿勢θsをθ1~θNと変化させて磁界を測定したときの、3次元磁気センサ60が順次出力する磁気データq1~qNを、センサ座標系ΣSにおいてプロットした図である。外部磁界Bxは不均一な磁界であり、地上座標系ΣGにおける位置に応じて、向き及び大きさが変化する。従って、3次元磁気センサ60の位置Psを変化させた場合には、3次元磁気センサ60が検知する外部磁界Bxの向き及び大きさも変化する。このような外部磁界Bxは、センサ座標系ΣSにおいて、3次元磁気センサ60の姿勢θs及び位置Psに依存して方向及び大きさの両方を変化させるベクトルSBxθs、Psとして表現される。そして、3次元磁気センサ60の姿勢θs及び位置Psを変化させつつ測定した複数の磁気データq1~qNは、内部磁界SBi、地磁気SBgθs、及び外部磁界SBxθs、Psのそれぞれのベクトルの和として出力される。従って、複数の磁気データq1~qNは、センサ座標系ΣSにおいて、中心点x0Gを中心とする球面SGと、球面とは異なる歪んだ形状を有する曲面SXとを、中心点x0Gを基準として重ね合わせた立体SDの表面近傍に分布することになる。ここで、曲面SXは、センサ空間ΣSにおいて外部磁界SBxθs、Psを表す曲面であり、3次元磁気センサ60が出力する複数の磁気データq1~qNの各々が示すベクトルのうち、内部磁界SBiと、外部磁界SBxθs、Psとの和の成分を表す。立体SDは、球面とは異なる形状を有するため、立体SDの表面近傍に分布する複数の磁気データq1~qNから、内部磁界SBiの値を示す中心点x0Gのベクトルを求めることは困難である。但し、外部磁界Bxが微弱である場合には、立体SDは球面に近い形状となる。図4Aは、外部磁界Bxが微弱である場合に、3次元磁気センサ60が順次出力する磁気データq1~qNを、センサ座標系ΣSにおいてプロットした図である。外部磁界Bxが微弱である場合、3次元磁気センサ60が順次出力する磁気データq1~qNは、地磁気を表すベクトルSBgθsの描く中心点x0Gを中心とする球面SGと、外部磁界を表すベクトルSBxθs、Psの描く曲面SXとの重ね合わせである、球面に近い形状の立体SDの表面近傍に分布する。この場合、複数の磁気データq1~qNから、中心点x0Gとほぼ一致するベクトルを求めることが可能である。なお、不均一な外部磁界Bxが存在する場合であっても、複数の磁気データq1~qNを取得する際の3次元磁気センサ60の位置PS及び姿勢θSの変化のさせ方によって、3次元磁気センサ60が、不均一な外部磁界Bxを、均一な磁界であるかのように検出する場合がある。具体的には、不均一な外部磁界Bxが存在する場合であっても、複数の磁気データq1~qNを取得する際に、利用者が機器1を手で握り、3次元磁気センサ60の位置PSが変化するように機器1を振るのではなく、機器1の位置厳密には、3次元磁気センサ60の位置PSを固定したままその姿勢のみを変化させる場合には、外部磁界Bxは、センサ空間ΣSにおいて、3次元磁気センサ60の姿勢θsに基づいて、その方向のみを変化させる一定の大きさのベクトルSBxθsとして表現される。この場合、3次元磁気センサ60から出力される複数の磁気データq1~qNは、図4Bに示すように、センサ空間ΣSにおいて、ほぼ球面の形状を有する立体SDの表面近傍に分布する。そして、立体SDの形状がほぼ球面となる場合、その中心点は、球面SGの中心点x0Gすなわち、内部磁界SBiを表す座標とほぼ等しくなる。この場合、複数の磁気データq1~qNから、中心点x0Gとほぼ一致するベクトルを求めることが可能である。このように、複数の磁気データq1~qNに基づいて、オフセットcすなわち、地磁気Bgを表す球面SGの中心点x0Gを算出するためには、立体SDが球面に近い形状である必要がある。つまり、外部磁界Bxの影響が大きく、立体SDの形状が球面とは大きく異なる歪んだ形状である場合には、複数の磁気データq1~qNに基づいて、オフセットcを算出することはできない。従って、外部磁界Bxの及ぼす影響の大きさ、すなわち、立体SDがどの程度球面に近い形状を有しているかを評価し、立体SDの形状がほぼ球面であると看做される場合のみ、当該球面立体SDの中心点をオフセットcとして採用し、立体SDの形状が球面とは異なる形状を有する場合には、複数の磁気データq1~qNに基づいた立体SDに基づいたオフセットcの算出を避ける必要がある。[2.磁気データが球面近傍に分布する場合の球面の中心点の算出]以下では、図5を参照しながら、3次元磁気センサ60が出力するN個の磁気データq1~qNが、半径rの球面中心点算出用球面Sの近傍に分布すると仮定する場合に、球面Sの中心点中心点算出用球面の中心点x0を算出する方法について述べる。ここで、球面Sは、センサ空間ΣS上で、複数の磁気データq1~qNが球面近傍に分布すると仮定して、その球面の中心点を算出するために便宜上導入される球面であり、前述した、地磁気を表す球面SG、及び立体SDとは異なるものである。なお、不均一な外部磁界Bxが微弱な場合、または、不均一な外部磁界Bxが存在する場合であっても機器1の位置厳密には、3次元磁気センサ60の位置を固定して姿勢のみを変化させつつ複数の磁気データq1~qNを測定した場合には、地磁気Bgを表す球面SG及び曲面SXを重ね合わせた立体SDと、球面Sとは、ほぼ一致する。具体的には、外部磁界Bxが微弱な場合、半径rは地磁気Bgの大きさを表し、球面S及びその中心点x0は、地磁気Bgを表す球面SG及びその中心点x0Gとほぼ等しくなる。図14C参照また、不均一な外部磁界Bxが存在する場合であっても、機器1の位置を固定して姿勢のみを変化させつつ複数の磁気データq1~qNを測定した場合、半径rは地磁気Bgを表すベクトルと外部磁界Bxを表すベクトルとを足し算したベクトルの大きさを表し、球面Sの中心点x0は、地磁気Bgを表す球面SGの中心点x0Gとほぼ等しくなる。図14B参照複数の磁気データq1~qNが、半径rの球面S上に存在すると仮定する場合、各磁気データqiと、球面Sの中心点x0との距離はrであるため、以下の式1~3が成立する。なお、以下では、特に断りの無い限り、ベクトルや座標等は、センサ空間ΣSにおいて表現されたものであるとする。ここで、複数の磁気データq1~qNの重心をqcとする。重心qcは、式4で表される。また、重心qcを原点とする重心座標系ΣCを導入する。センサ座標系ΣSにおいて表現された磁気データqi及び中心点x0と、重心座標系ΣCにおいて表現された磁気データCqi及び中心点Cx0との間には、以下の式6及び式7が成立する。センサ座標系ΣSにおける球面Sの方程式1は、重心座標系ΣCにおいて式8として表現される。そして、式8に対して、式6及び7を代入することで、式9が得られる。式9のqiに対して、複数の磁気データq1~qNの各々を代入した結果を、磁気データの個数Nで割り算することにより、以下の式10を得ることができる。そして、式10と、式9を展開した式11との差分を取ることで、変数qiに依存しない項を消去し、式12を得ることができる。このように、磁気データである変数qiが、球面S上に存在することを示す球面の公式である式1は、磁気データqi、重心qc、及び中心点x0から構成される式12に変形された。この式12の変数qiに対して、N個の磁気データq1~qNをそれぞれ代入することで得られるN個の方程式は、行列Xを用いて式13として表現される。式13は、全ての磁気データq1~qNの座標が、中心点x0を中心とする球面S上に完全に一致する場合には解を有する。しかし、3次元磁気センサ60の測定誤差等を考慮すると、全ての磁気データq1~qNの座標が、球面Sと完全に一致する位置に存在することは無いため、式13は解を持たない。そこで、統計的な手法により尤もらしい解を得るために、式17で表される誤差を吸収するベクトルである誤差中心点算出用球面誤差ベクトルδSを導入する。このような誤差δSのノルムを最小にするベクトルx、換言すれば、Tを最小化するようなベクトルxが、球面Sの中心点x0として尤もらしいものであるいえる。ここで、以下の式19で示される目的関数fSxを定義する。このとき、目的関数fSxを最小化するxが、球面Sの中心点x0として尤もらしい値となる。中心点x0は、式21において示される共分散行列Aが正則である場合には、式20により求めることができる。このように、複数の磁気データq1~qNが球面Sの近傍に確率的に分布すると仮定して、複数の磁気データq1~qNで示される3軸の座標と球面Sとの誤差を最小化することにより、球面Sの中心点x0を算出することができる。球面Sの中心点x0は、前述の通り、不均一な外部磁界Bxが微弱な場合、または、不均一な外部磁界Bxが存在する場合であっても機器1の位置を固定して姿勢のみを変化させつつ複数の磁気データq1~qNを測定した場合には、地磁気Bgを表す球面SGの中心点x0Gとほぼ一致する。すなわち、立体SDと、球面Sとがほぼ一致する場合には、球面Sの中心点x0は、地磁気を表す球面SGの中心点x0Gとほぼ一致し、球面Sの中心点x0を、オフセットcとして採用することが可能である。[3.磁気センサが検出する磁気データの分布]ところで、3次元磁気センサ60が出力する複数の磁気データq1~qNが球面Sの近傍に分布する場合に、球面Sの中心点x0を算出するためには、複数の磁気データq1~qNがセンサ座標系ΣSにおいて3次元的な広がりを有して分布していることが必要である。例えば図6に示すように、複数の磁気データq1~qNが、センサ座標系ΣSの平面π上の円πC近傍に2次元的に分布する場合、球面Sは、円πCを切断面に有するような球面であるということしか特定できない。円πCを切断面に有するような球面は、円πCの中心点πCOを通り平面πに直交する直線πL上に存在する中心点xπ1を中心とする球面Sπ1であるかもしれないし、直線πL上に存在する中心点xπ1とは異なる中心点xπ2を中心とする球面Sπ2であるかもしれない。つまり、中心点x0が、円πCの中心点πCOを通り平面πに直交する直線πL上に存在することまでは特定可能であっても、直線πL上のどの位置に存在するかについて、具体的に特定することが出来ない。機器1が自動車等のように左右方向に姿勢変化する場合や、機器1が携帯電話等のようにユーザの手によって動かされるがその動かし方が不十分な場合には、機器1の姿勢変化は3次元的とはならず、2次元的なものとなることがある。この場合、センサ座標系ΣSにおいて複数の磁気データq1~qNの各々が示す座標の分布は、3次元的な広がりを有さず、2次元的に広がるため、複数の磁気データq1~qNに基づいて中心点x0を算出することは出来ない。3次元磁気センサ60の測定誤差等の影響により、複数の磁気データq1~qNがわずかに3次元的な広がりを有する場合には、複数の磁気データq1~qNをその近傍に有するような何らかの球面Sπの中心点xπが算出される。しかし、この中心点xπは、複数の磁気データq1~qNの各々に含まれる測定誤差の影響によって確率的に決定された座標に過ぎない。中心点xπと、仮に、機器1の姿勢変化が3次元的であった場合に測定される複数の磁気データを近傍に有するような球面Sの中心点x0とは、異なる可能性が高い。例えば、図6の例のように、磁気データq1~qNが、平面πに垂直方向の誤差を有して分布する場合に、複数の磁気データq1~qNより球面Sπ1の中心点xπ1が算出されたとしても、実際には、球面Sπ2の中心点xπ2が、球面Sの中心点x0かもしれない。このような不正確な座標である中心点xπは、オフセットcとして採用することは出来ない。なお、図6のように、複数の磁気データq1~qNが2次元的に分布している場合であっても、機器1の姿勢変化が3次元的であった場合に測定される複数の磁気データを近傍に有するような球面Sの中心点x0が、直線πL上に存在することについては特定可能である。つまり、平面π上に2次元的に分布する複数の磁気データq1~qNに基づいて算出される中心点xπが有する情報のうち、平面πに垂直な方向つまり、直線πL方向の情報については正確性を有さないが、平面πに平行な方向の情報については正確性を有する。従って、オフセットcの算出にあたっては、中心点xπが有する情報のうち、正確性を有する方向の情報つまり、オフセットcの候補となる中心点x0が、直線πL上に存在しているという情報について、利用することが可能である。図7は、不均一な外部磁界Bxの影響が大きく、立体SDが球面とは大きく異なる形状を有する場合で、且つ、3次元磁気センサ60から順次出力された磁気データq1~qNが、平面π上に2次元的に分布する場合を示している。この場合、複数の磁気データq1~qNが、円周πC近傍に存在していると仮定すれば、球面Sの中心点x0は、円周πCの中心点πCOを通り平面πに直交する直線πL上に存在することを特定することができる。しかし、立体SDが球面とは大きく異なる歪んだ形状を有する場合には、立体SDの表面近傍に存在する複数の磁気データq1~qNは、球面Sの近傍には存在しない。従って、平面πによる球面Sの切断面である円周πC近傍に複数の磁気データq1~qNが存在していると仮定して算出された中心点x0と、地磁気Bgを表す球面SGの中心点x0Gとは一致しない可能性が高く、中心点x0が直線πL上に存在する場合であっても、中心点x0Gは直線πL上には存在しない可能性が高い。このように、複数の磁気データq1~qNが2次元的に分布し、球面Sの中心点x0が直線πL上に存在することを特定可能な場合であっても、外部磁界Bxの影響が大きく、地磁気Bgを表す球面SGの中心点x0Gと球面Sの中心点x0とが一致しない場合には、中心点x0を、オフセットcとして採用することを防止する必要がある。なお、複数の磁気データq1~qNが、平面π上に2次元的に分布する場合には、立体SDの3次元的な形状を知ることは困難である。しかし、立体SDが、球面とは大きく異なる歪んだ形状を有する場合、平面πによる立体SDの切断面を表す曲線πDも、円周とは異なる歪んだ形状を有する可能性が高い。つまり、複数の磁気データq1~qNが、立体SDの平面πによる切断面を表す曲線πDの近傍に分布していると仮定した場合、曲線πDの形状の歪の程度を評価することにより、立体SDの形状と球面の形状との相違の程度について判定することが可能である。そして、曲線πDの形状が円周とは異なる大きく歪んだ形状である場合には、立体SDと球面とが異なる形状を有し、地磁気Bgを表す球面SGの中心点x0Gと球面Sの中心点x0とが一致しない可能性が高い。従って、曲線πDが円周とは異なる歪んだ形状を有する場合には、球面Sの中心点x0が直線πL上に存在することが特定される場合であっても、球面Sの中心点x0を、オフセットcとして採用することを避ける必要がある。一方、曲線πDの形状が歪のない円周であると看做すことができる場合には、曲線πDと円周πCとがほぼ一致すると看做すことができる。この場合、立体SDと球面Sとがほぼ一致すると看做すことができ、球面Sの中心点x0を、オフセットcとして採用することができる。本実施形態では、複数の磁気データq1~qNが3次元的な広がりを持って分布している場合、及び複数の磁気データq1~qNが2次元的に分布している場合の双方について、複数の磁気データq1~qNに基づいて外部磁界Bxの影響の大きさ、すなわち、立体SDと球面との形状の相違の程度を評価し、当該評価結果に基づいて、球面Sの中心点x0をオフセットcに採用することの可否を判定する。以下、本実施形態における球面Sの中心点x0の算出方法、及び中心点x0のオフセットcへの採否の判定方法について、具体的に説明する。[4.機器構成及びソフトウェア構成]図8は、本発明の実施形態に係る機器1の構成を示すブロック図である。機器1は、各種の構成要素とバスを介して接続され装置全体を制御するCPU10、CPU10の作業領域として機能するRAM20、各種のプログラムやデータを記憶したROM30、通信を実行する通信部40、画像を表示する表示部50、及び地磁気を検出して地磁気データを出力する3次元磁気センサ60、を備える。3次元磁気センサ60は、X軸地磁気センサ61、Y軸地磁気センサ62、及びZ軸地磁気センサ63を備える。各センサは、MI素子磁気インピーダンス素子、MR素子磁気抵抗効果素子などを用いて構成することができる。地磁気センサI/F64は、各センサから出力信号をAD変換して磁気データを出力する。この磁気データは、X軸地磁気センサ61、Y軸地磁気センサ62、及びZ軸地磁気センサ63からの出力を、センサ座標系ΣSにおけるベクトルデータとして、センサ座標系ΣSの、x軸、y軸およびz軸の3成分によって示されるデータである。CPU10は、ROM30に格納されている磁気データ処理プログラム70を実行することによって、3次元磁気センサ60の出力からオフセットを算出すると共に、正確な地磁気の向き及び大きさを算出する。表示部50は、CPU10が磁気データ処理プログラム70を実行することにより算出した地磁気の向きを示す方位情報を、矢印等によって表示する。なお、磁気データ処理プログラム70は、地図アプリケーション等との連携を想定するものであっても良く、表示部50は、地磁気の向きを示す方位情報である矢印等を地図上に表示しても良い。CPU10、RAM20、3次元磁気センサ60、及び磁気データ処理プログラム70は、3次元磁気センサ60が検出する磁気データに基づき、正確な地磁気の向き及び大きさを示す地磁気データを算出する、地磁気測定装置として機能する。この地磁気データは、地上座標系ΣGにおいて磁極北に向かう水平成分と伏角方向の鉛直成分とを有するベクトルを、センサ座標系ΣSから見たベクトルデータとして、x軸、y軸およびz軸の3成分によって表されるデータである。磁気データ処理プログラム70は、バッファ管理モジュール71、オフセット算出モジュール72、及び方位算出モジュール73等のモジュール群で構成される。磁気データ処理プログラム70は、オフセットcを算出し、このオフセットc及び3次元磁気センサ60が検出する磁気データに基づいて、正確な地磁気データを算出する。なお、オフセットcは、常に一定の座標となるものではない。すなわち、オフセットcは、必要に応じて変更される必要がある。これは、機器1の内部状態が変化した場合、例えば、機器1に搭載された部品を流れる電流の大きさが変化した場合や、機器1に搭載された部品の着磁状況等の変化に伴い、内部磁界Biが変化するからである。このような内部磁界Biの変化を想定し、内部磁界Biを表すベクトルであるオフセットcを必要に応じて更新することによって、正確な地磁気Bgを測定することが可能となる。以下、説明の便宜上、オフセットcの更新にあたり、更新前のオフセットcを旧オフセットc0と呼び、更新後のオフセットcを新オフセットc1と呼ぶ。バッファ管理モジュール71は、3次元磁気センサ60から順次出力される磁気データq1~qNをバッファに蓄積するNは精度よくオフセットを導出するために必要な磁気データの規定測定回数を表す4以上の自然数。これら複数の磁気データq1~qNの蓄積手段としては、RAM20を用いる。オフセット算出モジュール72は、バッファ管理モジュール71に蓄積されている複数の磁気データq1~qNを用いて、新オフセットc1を算出し、オフセットcを、旧オフセットc0から新オフセットc1へと更新する。オフセットcの保持手段としては、RAM20を用いる。方位算出モジュール73は、3次元磁気センサから順次出力される磁気データq1~qNを、オフセット算出モジュール72によって保持されるオフセットcにより補正して、正確な地磁気データを算出し、方位情報を生成する。なお、方位算出モジュール73は、適当な周期でバッファ管理モジュール71及びオフセット算出モジュール72を呼び出し、オフセット更新の指示を与える。[5.処理の流れ]図9は、オフセットcを算出し更新する処理の流れを示すフローチャートである。このフローチャートは、方位算出モジュール73が、バッファ管理モジュール71及びオフセット算出モジュール72を呼び出しオフセット更新の指示を与えた場合に実行される。方位算出モジュール73は、上位プログラムからの指示がある場合等、オフセットcの更新が必要である場合に、オフセット更新の指示を行う。ステップS100において、CPU10は、初期化処理、及び磁気データ取得処理を行う。初期化処理において、CPU10は、バッファ管理モジュール71を呼び出し、RAM20に記憶した磁気データを廃棄する。なお、本実施形態では、磁気データの全部を廃棄するが、RAM20に蓄積されたデータのうち古い方から一定割合のデータのみを廃棄しても良い。そして、磁気データ取得処理において、3次元磁気センサ60から出力されるN個の磁気データq1~qNを、RAM20に格納する。ステップS101において、CPU10は、磁気データ分布指標算出処理を行う。磁気データ分布指標算出処理は、ステップS100でRAM20等に蓄積されたN個の磁気データq1~qNに基づいて、センサ座標系ΣSにおけるN個の磁気データq1~qNが示す座標の分布の広がりを示す各種分布指標を算出する処理である。具体的には、磁気データ分布指標算出処理では、N個の磁気データq1~qNの分散を表す共分散行列Aを算出すると共に、共分散行列Aの固有値λ1~λ3及び固有ベクトルu1~u3を算出する。これにより、N個の磁気データq1~qNの各々が示す座標が、センサ座標系ΣSにおける各固有ベクトルu1~u3方向のそれぞれについて、どの程度の広がりを有するか知ることができる。ステップS102において、CPU10は、中心点算出処理を行う。中心点算出処理は、ステップS100で取得したN個の磁気データq1~qNと、ステップS101で得た各種分布指標とに基づいて、N個の磁気データq1~qNを近傍に有する球面Sの中心点x0を算出する。中心点x0は、オフセットcを、旧オフセットc0から新オフセットc1に更新するときの、新オフセットc1の候補となる座標である。ステップS103において、CPU10は、歪判定処理を行う。歪判定処理は、ステップS100で取得したN個の磁気データq1~qN、ステップS101で算出した各種分布指標、及び、ステップS102で算出した中心点x0に基づいて、N個のデータの分布の形状が、前述した図2及び図4のような、ある点を中心とする球面の近傍に分布するものであるか否かについて判定する処理である。つまり、歪判定処理は、N個の磁気データq1~qNが、ある立体SDの表面近傍に分布すると仮定した場合に、立体SDの形状と、球面の形状との相違の程度を判定する処理である。そして、立体SDの形状が球面と看做せると判定した場合、すなわち、N個の磁気データq1~qNが図2及び図4に示すような球面に近い形状の立体の表面近傍に分布すると判定した場合には、CPU10は、処理をステップS104に進める。一方、立体SDが球面とは異なる大きく歪んだ形状を有すると判定した場合には、CPU10は、フローチャートに示される処理を終了する。ステップS104において、CPU10は、オフセット更新処理を行う。オフセット更新処理は、中心点x0を新オフセットc1として採用したうえで、オフセットcを、旧オフセットc0から新オフセットc1に更新する処理である。オフセットcを更新した後、CPU10は、フローチャートに示される処理を終了する。なお、本実施形態では、ステップS103において、立体SDが球面とは異なる大きく歪んだ形状を有すると判定した場合、CPU10は、フローチャートに示される処理を終了するが、本発明はこれに限定されるものではない。例えば、ステップS103において、立体SDが球面とは異なる大きく歪んだ形状を有すると判定した場合に、CPU10は、処理をステップS100に戻しても良いし、何らかのメッセージを表示部50に出力した上で処理を一旦停止させ、ユーザからの指示を待ってステップS100から処理を再開させても良い。N個の磁気データq1~qNを取得する際に、利用者が機器1を手で握って回転させるのではなく、機器1の位置を固定したままその姿勢のみを変化させるようにすると、外部磁界Bxの影響を低く抑えることができる図4B参照。そのため、ステップS103において、立体SDが球面とは異なる大きな歪を有する形状であると判定された場合には、機器1の位置を固定したまま回転させることを、利用者に対して指示しても良い。利用者に対する指示は、機器1の表示部50に画像や動画等を表示したり、音声等を用いたりすることで、行っても良い。以下において、これら、ステップS101~S104においてCPU10が行う、磁気データ分布指標算出処理、中心点算出処理、歪判定処理、及びオフセット更新処理の詳細について説明する。[6.磁気データ分布指標算出処理]図10を参照しつつ、ステップS101における、磁気データ分布指標算出処理について説明する。図10は、複数の磁気データq1~qNのセンサ座標系ΣSにおける、分布の様子を示す図であり、複数の磁気データq1~qNが3次元的に広がりを持って分布している場合を例示している。磁気データ分布指標算出処理では、式21に示される共分散行列Aを用いて、複数の磁気データq1~qNの分布が、どの程度3次元的な広がりを有しているかについての各種指標を算出する。以下で、共分散行列Aの性質を説明する。共分散行列Aの固有値を、大きい順番に最大固有値λ1、中間固有値λ2、最小固有値λ3とし、それぞれの固有値に対応する大きさ1に正規化された固有ベクトルをu1、u2、u3とする。また、磁気データqiを、前述した重心qCを原点とする重心座標系ΣCにおいて表したベクトルをCqiと表す。このとき、固有値λjj=1、2、3は、固有ベクトルuj方向の分散σ2jに等しい。図10に示すように、各固有ベクトルu1~u3を重心座標系ΣCの原点qCを起点となるように配置する。このとき、例えばj=1の場合について考える。固有値λ1は、ベクトルCqiを、固有ベクトルu1へ射影した長さLi1の二乗Li12を、N個の磁気データCqii=1、2、...Nについて平均した値に等しくなる。つまり、固有値λ1は、N個の磁気データCqiが、重心qCから固有ベクトルu1方向にどの程度離れているかを表す指標であり、複数の磁気データq1~qNの分布が固有ベクトルujの方向に対してどの程度の広がり有するかについて、対応する固有値λjの大きさに基づいて判断することができる。最小固有値λ3に対応する固有ベクトルu3の方向が、複数の磁気データq1~qNの分布の広がりが最も小さい方向であり、最小固有値λ3が、複数の磁気データq1~qNの分布の広がりが最も小さい方向における広がりの程度を示す指標である。従って、複数の磁気データq1~qNの分布の3次元的な広がりの程度は、最小固有値λ3の値によって評価できる。最小固有値λ3の値が大きければ、複数の磁気データq1~qNが3次元的な広がりを有して分布しているが、最小固有値λ3の値が小さい場合は、複数の磁気データq1~qNが2次元的に分布していることになる。このように、磁気データ分布指標算出処理では、式21に示される共分散行列Aの、固有ベクトルu3と、最小固有値λ3とを利用して、センサ座標系ΣSにおける複数の磁気データq1~qNの分布の広がりが最も小さな方向、及びその広がりの程度を算出する。[7.中心点算出処理]次に、ステップS102の中心点算出処理について説明する。前述の通り、内部磁界Biは機器1の内部状態が変化した場合に変化するものであるため、内部磁界Biを表すオフセットcも必要に応じて更新されることが必要となる。中心点算出処理は、このオフセットcを、旧オフセットc0から新オフセットc1に更新するときの、新オフセットc1の候補である中心点x0を算出する処理である。図11に示すように、ある期間において3次元磁気センサ60から出力される複数の磁気データq1~qNが平面π上に2次元的に分布している場合であっても、測定誤差の影響によりわずかに3次元的に分布していると看做される場合には、複数の磁気データq1~qNをその近傍に有するような何らかの球面Sπの中心点xπが算出される。このような中心点xπは、複数の磁気データq1~qNの各々が有する測定誤差の影響により確率的に算出された座標に過ぎないため、中心点xπと、当該期間において仮に機器1を3次元的に姿勢変化させていれば測定されたであろう複数の磁気データを近傍に有するような球面Sの中心点x0とは、異なるものである可能性が高い。しかし、球面Sの中心点x0が中心点xπを通り平面πに直交する直線πL上に存在することについては、特定可能である。つまり、中心点xπの座標が有する情報のうち、直線πLと平行な方向についての情報は不正確であるが、直線πLと垂直な方向についての情報は正確である。従って、当該期間が開始される前のオフセットc旧オフセットc0が直線πL上に存在しない場合には、オフセットcを旧オフセットc0から直線πL上の新オフセットc1に更新して補正処理を行うことで、オフセットcを旧オフセットc0から更新しない場合に比べて正確な地磁気Bgを算出することができる。中心点算出処理では、中心点xπで示される座標のうち正確性の高い方向の情報については中心点xπの情報を利用しつまり、当該期間において3次元磁気センサ60から出力される複数の磁気データq1~qNから算出される情報を利用し、中心点xπで示される座標のうち正確性の低い方向については旧オフセットc0の情報を利用することで、中心点x0を推定する。これにより、中心点xπの座標が有する情報のうち不正確な情報を排除し、正確な情報のみに基づいて、中心点x0を求めることができる。なお、旧オフセットc0から中心点x0を示すベクトルx0-c0を、更新ベクトル中心点算出用更新ベクトルkと呼ぶ。図11に示すように、複数の磁気データq1~qNが平面π上に2次元的に分布し、複数の磁気データq1~qNの分布が直線πL方向に広がりを有さない場合、直線πL方向と、共分散行列Aの固有ベクトルu3の示す方向とは一致する。従って、中心点xπの座標が有する情報のうち、固有ベクトルu3方向の情報は、不正確な情報である。また、平面πは、共分散行列Aの固有ベクトルu1及びu2の示す方向と平行である。従って、中心点xπの座標が有する情報のうち、固有ベクトルu1方向及びu2方向の情報は、正確な情報である。更新ベクトルkは、中心点xπの座標が有する情報のうち、正確性を有する方向の情報を利用して、中心点x0を算出するためのベクトルである。従って、複数の磁気データq1~qNが2次元的に分布する場合、更新ベクトルkは、固有ベクトルu3方向の成分を有さないように定められる。つまり、更新ベクトルkは、固有ベクトルu1及びu2の線形結合として表現される。このとき、更新ベクトルkと固有ベクトルu3との間には以下の式22が成立する。式22で示される制約条件は、複数の磁気データq1~qNの分布が、3次元的な広がりを有さずに完全に押しつぶされた形状となる場合、すなわち、固有ベクトルu3方向の広がりを全く有さず、完全に2次元的な分布となる場合を表す制約条件である。一方、複数の磁気データq1~qNの分布が、完全に押しつぶされること無く、固有ベクトルu3方向に広がりを有する場合、中心点xπの座標が有する情報は、固有ベクトルu3方向について完全に欠落している訳ではない。この場合、更新ベクトルkと固有ベクトルu3とが完全に直交する必要は無く、式22で示される制約条件は緩和されることが好ましい。そこで、本実施形態では、式22で示される制約条件を、磁気データq1~qNが示す座標の分布の、固有ベクトルu3方向の広がりの大きさつまり、共分散行列Aの最小固有値λ3の大きさに応じて緩和する。式22で示した制約条件を緩和するために、以下の式23に示す誤差第3誤差δUを導入する。誤差δUは、式24に示す内積ε3に対し、係数第1係数αを掛け算した値である。内積ε3は、更新ベクトルkと、固有ベクトルu3との内積である。内積ε3は、更新ベクトルkと、固有ベクトルu3とが直交する場合には0となるが、直交しない場合には0以外の値となる。なお、更新ベクトルkは、中心点x0の位置を求めるための変数であるベクトルxの示す位置を、旧オフセットc0を起点として表した位置ベクトルであるが、計算の便宜上、式13との整合を図り、式25に示すように、ベクトルx及び旧オフセットc0を重心qcから見た位置ベクトルとして表現する。式23に現れる係数αは、式22で示した制約条件を緩和するための係数である。誤差δUを、0または0に近い小さな値に近づけることで、更新ベクトルkを定めることを考える。更新ベクトルkと固有ベクトルu3とが直交しない場合には、内積ε3は大きな値を有することになるが、係数αを十分に小さな値例えば、0または0に近い小さな値に設定した場合、内積ε3が大きな値を有していても、誤差δUを0または0に近い小さな値に近づけることができる。一方、係数αを大きな値に設定した場合、更新ベクトルkと固有ベクトルu3とが直交し、内積ε3が0または0に近い小さな値にならない限り、誤差δUを0または0に近い小さな値に近づけることはできない。すなわち、式23において係数αを0または0に近い小さな値に設定して、誤差δUを0または0に近い小さな値に近づける場合、更新ベクトルkと固有ベクトルu3とが直交するという式22に示した制約条件は緩和される。一方、式23において係数αを大きな値に設定して、誤差δUを0または0に近い小さな値に近づける場合、式23は式22と同様の制約条件を与える式になる。なお、本実施形態において、係数αは、正の値を有し、以下の式26に示すような、固有値λ3の関数として定められる。なお、式26の定数kαは、kα>0を満たす実数である。係数αは、固有値λ3が大きな値の場合には、0に近い小さな値となり、固有値λ3が小さな値の場合には、1に近い大きな値となる。式17で示した誤差δS、及び23で示した誤差δUを最小化する中心点x0を算出するため、以下の式27に示す目的関数中心点算出関数gSxを導入する。目的関数gSxは、ベクトル中心点算出ベクトルδSUの大きさを表す関数である。ベクトルδSUは、式28に示すように、誤差δS及び誤差δUを要素とするN1次元のベクトルである。複数の磁気データq1~qNの分布が2次元的な場合、誤差δSを最小化することで、複数の磁気データq1~qNを近傍に有する球面Sπの中心点xπとして尤もらしい座標を算出できるが、これに加えて、誤差δUを最小化することで、中心点xπが示す座標の有する情報のうち正確性の高い方向の情報についてのみ利用した、中心点x0として尤もらしい座標を算出することができる。複数の磁気データq1~qNの分布が3次元的な場合、係数αが0または0に近い小さな値に設定されるため、誤差δUを最小化しても、式22に示した制約条件は緩和され、十分に機能しないまたは完全に機能しない。従って、目的関数gSxを最小化するxを算出することで、複数の磁気データq1~qNの分布が2次元的な場合及び3次元的な場合の双方の場合について、複数の磁気データq1~qNを近傍に有する球面Sの中心点x0として尤もらしい座標を求めることができる。以下において、目的関数gSxを最小化する解x=中心点x0を算出する方法について述べる。ベクトルδSUは、式28に対して、式17、及び式23~25を代入することで、式29のように変形される。ここで、行列XSは式30に示すような、N1行3列の行列であり、ベクトルjSは式31に示すN1次元のベクトルである。目的関数gSxを最小化するような解xとして求められる中心点x0は、式32に示す3元1次連立方程式の解として求められ、式35で示す3次元のベクトルとして表される。ここで、式32は、式27に対して、最小二乗法を適用することによって算出される。また、行列ASは式33に示す3行3列の行列であり、ベクトルXSTjSは、式34に示す3次元のベクトルである。[8.歪判定処理]次に、ステップS103の歪判定処理について説明する。なお、本実施形態の歪判定処理は、複数の磁気データq1~qNが3次元的に分布する場合、及び複数の磁気データq1~qNが2次元的に分布する場合の双方について、複数の磁気データq1~qNを表面近傍に有する立体SDの形状と球面の形状との相違の程度を判定する処理である。以下では、本実施形態に係る歪判定処理の前提として、複数の磁気データq1~qNが3次元的な広がりを有して分布する場合に限定した歪判定処理8.1節について説明したうえで、その後、本実施形態に係る歪判定処理について説明する8.2節。[8.1.磁気データが3次元的に分布する場合の歪判定処理]複数の磁気データq1~qNが3次元的に分布する場合で、且つ、不均一な外部磁界Bxの影響が小さく立体SDの形状がほぼ球面と看做せる場合、式20に示した球面Sの中心点x0は、地磁気を表す球面SGの中心点x0Gとほぼ一致する。しかし、図3に示したように、外部磁界Bxが一様ではなく、かつ、外部磁界Bxの影響が大きい場合には、3次元磁気センサ60が順次出力する複数の磁気データq1~qNは、中心点x0を中心とする球面Sの近傍には分布せず、球面とは異なる歪んだ形状を有する立体SDの表面近傍に分布することになる。この場合、球面Sの中心点x0は、地磁気Bgを表す球面SGの中心点x0Gとは異なる座標となる可能性が高いため、中心点x0をオフセットcとして採用することは、避ける必要がある図13参照。但し、立体SDと球面Sとが十分に近い形状である場合には、球面Sの中心点x0をオフセットcとして採用することができる。例えば、図4Aに示すように、外部磁界Bxの影響が小さい場合には、センサ座標系ΣSにおいて、3次元磁気センサ60が順次出力する複数の磁気データq1~qNが、ステップS102で算出された中心点x0を中心とする球面Sの近傍に分布することになり、且つ、当該球面S及びその中心点x0は、地磁気Bgを表す球面SG及びその中心点x0Gとほぼ等しくなる図14C参照。また、図4Bに示すように、機器1の位置厳密には、3次元磁気センサ60の位置を固定したままその姿勢のみを変化させて複数の磁気データq1~qNを取得する等、外部磁界Bxが3次元磁気センサ60によってほぼ一様な磁界として検出される場合には、複数の磁気データq1~qNが、ステップS102で算出された中心点x0を中心とする球面Sの近傍に分布することになる。そして、この場合、球面Sは地磁気Bgを表す球面SGとは異なるが、球面Sの中心点x0は、球面SGの中心点x0Gとほぼ等しくなる図14B参照。本節に係る歪判定処理は、複数の磁気データq1~qNが立体SDの表面近傍に分布すると仮定した場合に、立体SDの形状が球面の形状とどの程度近いものであるかを評価する処理である。立体SDの形状が球面の形状に十分に近いと評価された場合には、球面Sの中心点x0を、新オフセットc1として採用することができる。前述した中心点算出処理では、複数の磁気データq1~qNが、式13で示される球面Sの近傍に分布すると仮定した場合における、複数の磁気データq1~qNの各々が示す座標と球面Sとの複数の誤差を、誤差δSを用いて表した。これに対して本節に係る歪判定処理では、図12に示すように、まず、複数の磁気データq1~qNが、球面とは異なる歪んだ形状の式36に示される立体SDの表面近傍に分布すると仮定する。式36に示される立体SDは、球面を表す成分Xx-qc-jと、式37に示す歪を表す成分歪誤差ベクトルkEとを、加算したものである。ここで、式36の右辺に現れる0Nの右下の添え字Nは、零ベクトルの次元を表している。なお、式36に現れる成分Xx-qc-jが表す球面と、中心点算出処理により算出される球面Sとは、必ずしも一致しない。そこで、以下では、式36に現れる成分Xx-qc-jによって表される球面を球面歪判定用球面S2と称し、その中心点を中心点歪判定用球面の中心点x02と称する。なお、本節において、ベクトルxは、球面S2の中心点x02を表すための変数である。歪判定処理は、式36のうち、歪を表す成分kEの大きさを評価することにより、立体SDの形状と、球面S2の形状とが、どの程度相違しているかについて評価するものである。式37に示した歪を表す成分kEはN次元のベクトルである。歪を表す成分kEの第i行目の要素keは、以下の式39で表される関数kepに対して、球面Sの中心点x0を原点とする座標系から見た磁気データqiを表すベクトルqi-x0を代入することにより与えられる値である。関数kepは、式38に示す歪評価行列Eを係数行列とし、式40に示される3変数p=px、py、pzの二次形式で表される関数である。なお、歪評価行列Eは、式38に示したように、3行3列の対称行列である。なお、3次元磁気センサ60の測定誤差等を考慮すると、複数の磁気データq1~qNの各々が示す座標の全てが、立体SDと完全に一致する位置に存在することは無いため、式36は解を持たない。そこで、統計的な手法により尤もらしい解を得るために、式41で表される誤差を吸収するベクトルである誤差立体誤差ベクトルδSDを導入する。誤差δSDは、複数の磁気データq1~qNの各々が示す座標と立体SDの表面との誤差を各要素とする、N次元のベクトルであり、誤差歪判定用球面誤差ベクトルδS2と歪を表す成分kEとの和として表現される。なお、詳細は後述するが、誤差δS2は、式36のうち、球面を表す成分Xx-qc-jに対応する成分である。誤差δSDのノルムを最小にするベクトルx及び歪評価行列E、換言すれば、Tを最小化するようなベクトルx及び歪評価行列Eによって、複数の磁気データq1~qNを表面近傍に有する立体SDが表現される。ここで、以下の式42で示される目的関数立体歪評価関数fSDE,xを定義する。この目的関数fSDE,xを最小化するようなx及び歪評価行列Eにより、複数の磁気データq1~qNを近傍に有する立体SDが表現される。歪を表す成分kEの大きさすなわち、立体SDの形状と、球面S2の形状との、相違の程度については、後述する式43及び44で示される歪評価値gDEに基づいて評価することができる。そして、歪評価値gDEが閾値歪許容値δ0以下であれば、立体SDと球面S2とは等しい図形であると看做すことができるため、立体SDの形状は球面Sともほぼ等しくなり、中心点算出処理で算出した中心点x0をオフセットcとして採用することができる。一方、歪評価値gDEが閾値δ0よりも大きな値を有する場合には、中心点x0をオフセットcとして採用することはできない。以下で、式17で表される誤差δSと、式41で表される誤差δSDとの、性質の相違点について説明する。式17に表される誤差δSは、N個の磁気データq1~qNの各々が示す座標と球面Sとの誤差を吸収するためのベクトルである。誤差δSを構成する1行目~N行目の各要素は、それぞれ独立した変数である。従って、誤差δSによって、N個の磁気データq1~qNの各々が示す座標と球面Sとの誤差を吸収する場合には、N個の磁気データq1~qNと球面SとのN個の誤差の各々は、互いに制約の無い独立した値となる。つまり、誤差δSによって表されるN個の誤差は、それぞれが独立に確率的に定められるものであり、N個の誤差は全体として、対称性を有し、かつ方向依存性の無いホワイトノイズである。従って、式19乃至式21は、N個の磁気データq1~qNの座標と球面Sとの誤差を、対称性を持ち方向依存性の無いホワイトノイズとみなし、そのホワイトノイズである誤差を最小化するような球面Sの中心点x0を求めるものである。一方、式41に表される誤差δSDは、N個の磁気データq1~qNの各々が示す座標と立体SDとの誤差を吸収するベクトルであり、換言すれば、N個の磁気データq1~qNの各々が示す座標と球面S2との誤差を、誤差δS2と、歪を表す成分kEと、によって吸収するベクトルである。誤差δS2は、誤差δSと同様に式17の右辺で表され、N個の磁気データq1~qNの各々が示す座標と、球面S2との誤差を、ホワイトノイズとして表現するベクトルである。一方、歪を表す成分kEは、式39で示した3変数二次形式の関数kepを各成分とするベクトルである。3変数の二次形式は、変数が2次の項から構成される関数であり、3次元空間上の様々な曲面、例えば、直線、平面、柱面、球面、楕円面、錐面、1葉双曲面、2葉双曲面、及び各種放物面等を描くことができる。従って、歪を表す成分kEは、N個の磁気データq1~qNが示すN個の座標と球面S2とのN個の誤差の各々を、互いに独立したものではなく、N個の誤差の全てが同一の関数kepにより表される3次元空間内の曲面上に存在するという制約を持った値として表現する。このように、誤差δSDは、N個の磁気データq1~qNの各々が示す座標と球面S2とのN個の誤差の全体を、ホワイトノイズである誤差δS2と、滑らかな歪を表す成分kEによる球面S2からの歪という性質を有する誤差とによって表現する。つまり、誤差δSDを導入することにより、N個の磁気データq1~qNの各々が示す座標と球面S2とのN個の誤差を表すベクトルを、ホワイトノイズに起因する誤差δS2と、歪を表す成分kEによる歪に起因する誤差とに分離することが可能となった。以下では、球面Sの中心点x0を新オフセットc1として採用することが可能な場合と、中心点x0を新オフセットc1として採用することが不可能な場合とが、歪を表す成分kEの大きさに応じてどのように場合分けされるかについて、図13及び図14を参照しながら説明する。図13は、歪を表す成分kEが大きい場合、すなわち、不均一な外部磁界Bxの影響により、立体SDが球面S2とは大きく異なる形状を有する場合を示す。なお、図13及び図14では、簡単のため、複数の磁気データq1~qNが立体SDの表面に完全に一致するように分布すると仮定し、式36が成立する場合を示している。図13に示す場合、複数の磁気データq1~qNは、球面S2とは大きく異なる歪んだ形状の立体SD上に分布するため、複数の磁気データq1~qNが球面S上に分布していると看做すこともできない。中心点算出処理は、球面Sの近傍に複数の磁気データq1~qNが存在していることを前提として、オフセットc地磁気を表す球面SGの中心点x0Gの候補である中心点x0を算出する処理であるため、複数の磁気データq1~qNを近傍に有さない球面Sの中心点x0は、オフセットcとしての意味を有さない。また、立体SDは、地磁気Bgを表す球面SGと外部磁界Bxを表す曲面SXとの重ね合わせであると同時に、球面S2と歪を表す成分kEとの重ね合わせでもある。歪を表す成分kEは、複数の磁気データq1~qNの各々が示す座標と球面S2との誤差のうちホワイトノイズとして表現できない誤差を表すベクトルであり、外部磁界Bxを表す曲面SXとは異なるものである。よって、歪を表す成分kEが大きい場合には、球面S2と球面SGとは異なるものとなる可能性が高く、仮に球面Sと球面S2とが一致すると看做すことができる場合であっても、球面S及びその中心点x0と球面SG及びその中心点x0Gとは異なるものとなる可能性が高い。このように、式36が成立し、且つ、歪を表す成分kEが大きい場合には、球面Sの中心点x0と球面SGの中心点x0Gとが一致しない可能性が高く、球面Sの中心点x0をオフセットcとして採用することを避ける必要がある。一方、歪を表す成分kEが0となる場合には、球面Sの中心点x0をオフセットcとして採用することができる。例えば、図14Aに示すような、外部磁界Bxが存在しない場合を検討する。この例では、地磁気Bgを表す球面SGと外部磁界Bxを表す曲面SXとの重ね合わせである立体SDは、球面SGと等しくなり、その形状は球面となる。立体SDの形状が球面となり、歪を表す成分kEが0となる場合、立体SDと球面S2とは一致するため、式36が成立する場合には、式13も同時に成立する。つまり、球面Sは、複数の磁気データq1~qNを近傍に有するように定められる球面であるため、立体SDと球面Sとは一致する。その結果、球面Sと、球面SGとが一致することになるため、球面Sの中心点x0は、中心点x0Gを表すものとなり、中心点x0をオフセットcとして採用することができる。また、歪を表す成分kEが小さい場合歪を表す成分kE≒0となる場合には、球面Sの中心点x0をオフセットcと看做して採用することができる。例えば、図14Bに示すように、外部磁界Bxが3次元磁気センサ60によってほぼ一様な磁界として検出される場合を検討する。この例では、外部磁界Bxを表す曲面SXの形状を、中心点x0Gを中心とする球面と看做すことができるため、球面SGと曲面SXとの重ね合わせである立体SDは、中心点x0Gを中心とする球面に近い形状を有する。立体SDが球面に近い形状となり、歪を表す成分kEが0に近い小さな値となる場合、立体SDと球面S2とは、ほぼ等しくなるため、式36が成立する場合には、式13も成立すると看做すことができる。よって、立体SDと球面Sとは、ほぼ等しいものとなり、立体SD上に分布する複数の磁気データq1~qNは、中心点x0を中心とする球面S上にも分布していると看做すことができる。その結果、球面Sの中心点x0と、球面SGの中心点x0Gとは、ほぼ同一の座標になる。従って、図14Bに示すような場合にも、球面Sの中心点x0をオフセットcとして採用することができる。また、図14Cに示すような、3次元磁気センサ60によって検出される不均一な外部磁界Bxの大きさが十分に小さく、曲面SXの影響を無視することができる場合について検討する。この例では、地磁気Bgを表す球面SGと外部磁界Bxを表す曲面SXとの重ね合わせである立体SDは、地磁気Bgを表す球面SGとほぼ等しいと看做すことができるため、立体SDの形状は球面と看做すことができる。立体SDが球面に近い形状となり、歪を表す成分kEが0に近い小さな値となる場合、立体SDと球面S2とは、ほぼ等しくなるため、式36が成立する場合には、式13も成立すると看做すことができ、立体SDと球面Sとがほぼ等しいものとなる。よって、球面Sと球面SGとは、ほぼ等しくなるため、中心点x0と、球面SGの中心点x0Gとは、ほぼ同一の座標になる。従って、図14Cに示すような場合にも、球面Sの中心点x0をオフセットcとして採用することができる。このように、歪を表す成分kEが小さい場合、立体SD上に分布する複数の磁気データq1~qNは、球面S上にも分布すると看做すことができ、球面Sの中心点x0をオフセットcとして採用することができる。以下では、歪評価行列Eの性質について述べたうえで、歪を表す成分kEの大きさを評価する方法について述べる。前述の通り、N次元のベクトルである歪を表す成分kEの各要素は、3行3列の対称行列である歪評価行列Eを係数とする3変数二次形式で表される。すなわち、歪を表す成分kEのうち、第i行目の要素keqi-x0は、球面Sの中心点x0を起点として磁気データqiの座標を表したベクトルqi-x0と、ベクトルqi-x0を歪評価行列Eにより変換して得られるベクトルEqi-x0との内積である。ここで、歪評価行列Eの固有値を、絶対値の大きい順番に最大固有値λE1、中間固有値λE2、最小固有値λE3とし、それぞれの固有値に対応する大きさ1に正規化された固有ベクトルをuE1、uE2、uE3とする。歪評価行列Eは対称行列であるため、歪評価行列Eの3つの固有ベクトルuE1、uE2、uE3は互いに直交する。ベクトルqi-x0が、固有ベクトルuE1と平行な場合、ベクトルEqi-x0は、ベクトルqi-x0に平行で、ベクトルqi-x0をλE1倍した長さを有するベクトルとなる。つまり、要素keqi-x0の値は、λE1|qi-x0|2となる。従って、要素keqi-x0の絶対値は、ベクトルqi-x0の表す方向が、固有ベクトルuE1の表す方向に近づくに従って大きくなる。歪評価行列Eは、ベクトルxと歪評価行列Eとを変数とする目的関数fSDE,xを最小化することで求められる行列である。なお、前述の通り、本節においてベクトルxは、球面S2の中心点x02を表すための変数である。ここで、立体SDの表面近傍に存在する磁気データqiの座標が、球面S2から大きく離れた位置である場合、磁気データqiの示す座標と球面S2との誤差を、球面S2からのホワイトノイズである誤差δS2と看做すことは好ましくない。歪を表す成分kEは、立体SDの近傍に存在する複数の磁気データq1~qNの各々が示す座標と球面S2との誤差のうち、ホワイトノイズと看做すことが不適切な大きな誤差を、歪に起因する誤差として吸収するために定められるベクトルである。つまり、歪を表す成分kEは、複数の磁気データq1~qNを表面近傍に有する立体SDと球面S2との形状の相違を表すように定められる。そして、歪評価行列Eの各成分は、磁気データqiの示す座標が球面S2から大きく離れている場合には、歪を表す成分kEのうち第i行目の要素keqi-x0が大きな値となるように定められる。具体的には、歪評価行列Eの最大固有値λE1に対応する固有ベクトルuE1の方向が、ベクトルqi-x0の方向に近づくように、歪評価行列Eの各成分が定められる。なお、歪を表す成分kEは、複数の磁気データq1~qNと球面S2とのN個の誤差を全体として表すベクトルである。従って、歪評価行列Eの各成分は、歪評価行列Eの最大固有値λE1に対応する固有ベクトルuE1の方向と、球面S2からの誤差が大きな磁気データqiが多く存在する位置を球面Sの中心点x0から見た方向とが近づくように、定められる。そして、歪評価行列Eの最大固有値λE1は、中心点x0から見て、固有ベクトルuE1の方向に存在する複数の磁気データと球面S2との誤差の大きさの程度を表す値となる。このように、歪評価行列Eによって、立体SDと球面S2との形状の相違の程度と、形状の相違が大きな方向とを表現することが可能となる。ここで、以下の式43及び44に示す歪評価値gDEを、歪を表す成分kEの大きさを評価する評価値として定義する。歪評価値gDEは、目的関数fSDE,xを最小化する場合の、歪評価行列Eのノルム、すなわち歪評価行列Eの最大固有値λE1の絶対値である。歪評価値gDEの値が0となる場合には、式36は式13と一致し、立体SDが球面Sと等しくなる。また、歪評価値gDEの値が0とならない場合であっても、一定の閾値δ0以下であれば、立体SDが球面Sに近い形状を有するものであると判断できる。これらの場合には、複数の磁気データq1~qNが球面Sの近傍に分布すると看做せるため、前述した中心点算出処理により求められた球面Sの中心点x0をオフセットcとして採用することができる。以下に、歪評価値gDEを求める方法について説明する。まず、式39に示した関数kepは、二次形式で表される関数であるため、以下の式45のように、ベクトルの内積として表現することもできる。また、歪を表す成分kEの第i行目の要素keは、式39の変数pに対して、ベクトルを代入した値であるため、式47によって示される6次元のベクトルke2i、及び、式48によって表される評価行列Eの各成分を並べた6次元のベクトルeEを用いて、式46のように表せる。ここで、式49で示す行列X2を導入する。行列X2は、第i行目に、式47によって示される6次元のベクトルke2iを転置したものと、3次元のベクトルを転置したものを並べた、N行9列の行列である。式49で表される行列X2を用いることで、式42に示される目的関数fSDE,xを、式50で表される目的関数gSDeに変形することができる。なお、ベクトルeは、式51に示すように、式48に示した6次元のベクトルeEと、式52に示した3次元のベクトルeXとを並べた、9次元のベクトルである。また、式52のベクトルeXは、球面S2の中心点x02を表すための変数であるベクトルxから、式4に示した重心qcを引いたベクトルである。式50で表される目的関数gSDeを最小化する解e=e0は、以下の式53式で表現される連立方程式に対して、ガウス消去法や、コレスキー分解法などを適用することで、求められる。なお、式53は、式50に対して、最小二乗法を適用することによって算出される。このようにして得た解e0に基づき、式38により歪評価行列Eを復元する。そして、式43に示した歪評価値gDEの値、即ち、歪評価行列Eのノルムを求める。歪評価行列Eのノルムは、歪評価行列Eが有する3つの固有値のうち絶対値が最大となる固有値最大固有値λE1の絶対値に等しいため、ヤコビ法や、冪乗法により求めることができる。このとき、球面S2の中心点x02は、目的関数gSDeを最小化する解e0のうち式51のeXに相当する部分を式52に代入することにより求められる。歪評価値gDEの値が閾値δ0以下である場合には、中心点算出処理により求められた球面Sの中心点x0と、球面S2の中心点x02とは、ほぼ等しい座標となる。従って、歪評価値gDEの値が閾値δ0以下である場合には、球面Sの中心点x0、または、球面S2の中心点x02のいずれについても、オフセットcとして採用することができる。以上で述べた、本節に係る歪判定処理は、複数の磁気データq1~qNの分布が3次元的な広がりを有していることを前提とする。しかし、複数の磁気データq1~qNが2次元的に分布する場合には、式49からも明らかなように、式53に示される9行9列の行列X2TX2の階数が、9よりも小さくなる場合がある。この場合には、式53に基づいて、解e0を算出することはできず、評価行列Eを求めることもできない。本実施形態は、複数の磁気データq1~qNが3次元的に分布する場合のみならず、複数の磁気データq1~qNが2次元的に分布する場合も想定する。従って、複数の磁気データq1~qNが2次元的に分布する場合にも、複数の磁気データq1~qNを表面近傍に有するような立体SDが、球面からどの程度歪んだ形状を有するものであるかについての評価を行うことが必要である。そこで、以下、8.2節において、複数の磁気データq1~qNが2次元的に分布する場合にも、立体SDの形状と球面の形状との相違の程度を評価する方法について述べる。[8.2.本実施形態に係る歪判定処理]以下で、本実施形態における歪判定処理について述べる。本実施形態に係る歪判定処理は、8.1節で述べた、複数の磁気データq1~qNが3次元的に分布する場合にのみ適用可能な歪判定処理を、複数の磁気データq1~qNが2次元的に分布する場合にも適用できるように拡張するものである。図7のように、複数の磁気データq1~qNが平面π上に2次元的に分布する場合には、複数の磁気データq1~qNを表面近傍に有するような立体SDの3次元的な形状を正確に表現することはできない。しかし、この場合であっても、立体SDの平面πによる切断面である曲線πDの2次元的な形状を表現することはできる。従って、複数の磁気データq1~qNを近傍に有するような平面π上の曲線πDの形状と、平面π上の円周の形状との相違の程度を評価することで、立体SDと球面との形状の相違の程度を評価することができる。前述の通り、複数の磁気データq1~qNが平面π上に2次元的に分布する場合には、複数の磁気データq1~qNの分散を表す共分散行列Aの固有ベクトルu1及びu2は、平面πに平行になり、固有ベクトルu3は、平面πに垂直になる。すなわち、平面π上の曲線πDと、平面π上の円周との形状の相違の程度を評価するためには、固有ベクトルu1及びu2方向についての、立体SDと球面S2との形状の相違の程度を評価すればよい。なお、複数の磁気データq1~qNが、平面π上に2次元的に分布する場合には、固有ベクトルu3方向について、立体SDの形状を知ることはできない。従って、固有ベクトルu3方向について、立体SDと球面S2との形状の相違の程度を評価することは避けることが望ましい。そこで、本実施形態では、共分散行列Aの最小固有値λ3に対応する固有ベクトルu3方向については、立体SDと球面S2との形状の相違が存在しないと仮定して、歪評価行列Eを定める。具体的には、歪評価行列Eの最小固有値λE3が0となり、且つ、共分散行列Aの最小固有値λ3に対応する固有ベクトルu3と歪評価行列Eの最小固有値λE3に対応する固有ベクトルuE3とが平行となるという制約の下で、歪評価行列Eを定める。このような歪評価行列Eを定めた場合、共分散行列Aの固有ベクトルu3を歪評価行列Eによって変換したベクトルEu3は0となり、共分散行列Aの固有ベクトルu3とベクトルEu3との内積も0となる。よって、固有ベクトルu3とベクトルqi-x0とが平行となるような磁気データqiについて、歪を表す成分kEの要素keが、0となる。すなわち、本実施形態に係る歪評価行列Eは、中心点x0から見て共分散行列Aの固有ベクトルu3方向に存在する磁気データqiと球面S2との誤差が0であると看做す。換言すれば、歪評価行列Eは、固有ベクトルu3方向について、立体SDと球面S2との形状の相違が存在しないことと看做す行列として定められる。以上に述べたような、共分散行列Aの固有ベクトルu1及びu2方向については立体SDと球面S2との形状の相違を評価するが、固有ベクトルu3方向については立体SDと球面S2との形状の相違が存在しないことと看做す歪評価行列Eを算出するための制約条件第2の制約条件を、式54に示す。式54は、共分散行列Aの固有ベクトルu3を歪評価行列Eによって変換したベクトルEu3が、零ベクトルであることを表す。歪評価行列Eは対称行列であるため、歪評価行列Eの3つの固有ベクトルuE1~uE3は、互いに直交する。そして、共分散行列Aの固有ベクトルu3と、ベクトルEu3とは、共に、歪評価行列Eの固有ベクトルuE1~uE3の線形結合として表される。例えば、u3=v1uE1v2uE2v3uE3と表し、Eu3=w1uE1w2uE2w3uE3と表すことができる。ある行列の固有ベクトルに平行なベクトルを、当該ある行列により変換した場合、変換後のベクトルは、変換前のベクトルと平行で、且つ、当該ある行列の固有ベクトルに対応する固有値倍の大きさを有する。よって、固有ベクトルu3と、ベクトルEu3との間には、w1=λE1v1、w2=λE2v2、及びw3=λE3v3という関係が成立する。従って、式54は、歪評価行列Eの3つの固有値λE1~λE3の全てが0となること、または、共分散行列Aの固有ベクトルu3が、歪評価行列Eの3つの固有値λE1~λE3のうち値が0となる固有値に対応する固有ベクトルuE1~uE3の線形結合によって表されることを示す。複数の磁気データq1~qNが平面π上に2次元的に分布し、且つ、平面π上において立体SDと球面S2との形状の相違が存在する場合は、歪評価行列Eの最大固有値λE1及び中間固有値λE2は、0よりも大きな値となる。従って、式54は、歪評価行列Eの最小固有値λE3が0であり、且つ、共分散行列Aの固有ベクトルu3と歪評価行列Eの固有ベクトルuE3とが平行であることを表すことになる。なお、平面π上において立体SDと球面S2との形状の相違が存在しない場合、つまり、図7において、曲線πDと円周πCとが一致する場合、式54は、歪評価行列Eの3つの固有値λE1~λE3が全て0となることを表す。このように、式54に表現される制約条件を導入することにより、歪評価行列Eは、固有ベクトルuE1及びuE2が平面πと平行になり、固有ベクトルuE3が直線πLと平行になるように定められる。そして、歪評価行列Eは、共分散行列Aの固有ベクトルu1及びu2方向についてのみ、立体SDと球面S2との形状の相違の程度を評価し、固有ベクトルu3方向については、歪が無いものと看做す行列としての性質を有することになる。なお、本実施形態では、複数の磁気データq1~qNが2次元的に分布している場合のみならず、8.1節で述べたような、複数の磁気データq1~qNが3次元的な広がりを有して分布している場合であっても、立体SDと球面S2との形状の相違の程度を評価する必要がある。そして、複数の磁気データq1~qNの分布が3次元的な広がりを有する場合には、共分散行列Aの固有ベクトルu3方向についても、立体SDと球面S2との形状の相違の程度を評価する必要がある。従って、立体SDと球面S2との形状の相違の程度を評価するにあたっては、式54に示される制約条件を常に厳格に満たすべきではなく、式54に示される制約条件を、複数の磁気データq1~qNの3次元的な広がりの程度に応じて緩和する必要がある。そこで、式55に示すような、誤差第2誤差δEを導入する。誤差δEは、式54の左辺に表されるベクトルEu3に対して、係数第2係数βを掛け算した、3次元のベクトルである。係数βは、複数の磁気データq1~qNが2次元的に分布する場合には、大きな値となり、複数の磁気データq1~qNが3次元的に分布する場合には、0または0に近い小さな値を有する。ここで、誤差δEを最小化して0または0に近い小さな値に近づける場合を想定する。係数βが大きな値を有する場合、式55のうち、ベクトルEu3自体が0に近づかない限りは、誤差δEを0に近づけることはできない。すなわち、係数βが大きな値を有する場合には、式55は、式54に示される制約条件と等しい条件を示す式となる。一方、係数βが0または0に近い小さな値を有する場合、ベクトルEu3自体が大きな値を有していても、誤差δEを0に近づけることができる。すなわち、係数βが0または0に近い小さな値を有する場合には、式55は、式54に示される制約条件を緩和した条件を示す式となる。このように、誤差δEを最小化して0に近づけることで、複数の磁気データq1~qNの分布の3次元的な広がりの程度に応じた、立体SDの形状と球面S2の形状との相違の程度を評価することが可能となる。本実施形態において、係数βは、以下の式56に示すような、正の値を有する、固有値λ3の関数として定められる。なお、式56の定数kβは、kβ>0を満たす実数である。係数βは、固有値λ3が大きな値の場合には、0に近い小さな値となり、固有値λ3が小さな値の場合には、1に近い大きな値となる。式55に示した誤差δEは、式58に表される行列Mによって、式57のように変形される。ここで、行列Mは、式59に示す共分散行列Aの最小固有値λ3に対応する固有ベクトルu3の各要素を成分に有する、6行3列の行列である。また、式57の03×3は、3行3列の零行列である。ところで、8.1節磁気データが3次元的に分布する場合の歪判定処理に示した歪判定処理は、目的関数fSDE,xを最小化することで、歪評価行列Eの各成分と、球面S2の中心点x02とを算出し、歪評価行列Eのノルムに基づいて、複数の磁気データq1~qNを近傍に有する立体SDと、球面S2との形状の相違の程度を評価するものであった。目的関数fSDE,xを最小化することにより球面S2の中心点x02を算出する方法は、歪を表す成分kEを考慮しつつ球面を表す成分Xx-qc-jを最小化するものであり、2節磁気データが球面近傍に分布する場合の球面の中心点の算出に示した方法、つまり、式17に示した誤差δSを最小化することにより球面Sの中心点x0を算出する方法を利用するものである。従って、球面S2の中心点x02は、複数の磁気データq1~qNが3次元的に分布することを前提として算出される。しかし、複数の磁気データq1~qNが2次元的に分布する場合には、複数の磁気データq1~qNが3次元的に分布することを前提とする2節に示した方法に基づいて球面S2の中心点x02を算出することはできない。そこで、複数の磁気データq1~qNが2次元的に分布する場合の歪判定処理では、7節中心点算出処理に示した、複数の磁気データq1~qNが2次元的に分布する場合に球面の中心点を算出するための制約条件である式22を考慮して、歪判定処理で用いる球面S2の中心点x02を算出する。具体的には、歪評価行列Eの各成分と、球面S2の中心点x02とを算出する際に、式23に示された誤差δUを最小化するという条件を付与することで、複数の磁気データq1~qNの3次元的な広がりの程度に応じて、球面S2の中心点x02を算出する。なお、式22に示した更新ベクトルkは、球面Sの中心点x0を算出するためのベクトルであるが、歪判定処理では、球面S2の中心点x02を算出する必要がある。そこで、旧オフセットc0から球面S2の中心点x02を示すベクトルを更新ベクトル歪判定用更新ベクトルk2と称し、更新ベクトルk2と固有ベクトルu3とが直交する条件第1の制約条件を緩和するための誤差を、誤差第1誤差δU2と称する。式23に基づいて誤差δU2は、以下の式60のように表される。ここで、式61に示すベクトル歪評価ベクトルδSD2を導入する。ベクトルδSD2は、式41に示した誤差δSD、式55に示した誤差δE、及び式60に示した誤差δU2を要素に持つ、N4次元のベクトルである。なお、式41に示した誤差δSDは、計算の便宜上、式50に現れるベクトルX2e-jと変形したうえで、式61に適用する。同様に、式55に示した誤差δEは、計算の便宜上、式57に示した式に変形したうえで、式61に適用する。ベクトルδSD2を最小化する場合、誤差δSD、誤差δU2、及び誤差δEが全体として最小化される。誤差δU2を最小化することで、複数の磁気データq1~qNの分布が2次元的な場合及び3次元的な場合の双方の場合において、球面S2の中心点x02を算出することが可能となる。これにより、複数の磁気データq1~qNの分布が2次元的な場合であっても、歪評価行列Eを算出することができる。また、誤差δEを最小化することで、複数の磁気データq1~qNの分布が広がりを有する方向についてのみ立体SDと球面S2との形状の相違の程度を評価する歪評価行列Eを算出することが可能となる。このように、誤差δU2、及び誤差δEを含むベクトルδSD2を最小化することにより算出される歪評価行列Eによって、複数の磁気データq1~qNの分布の広がりが2次元的な場合及び3次元的な場合の双方の場合において、立体SDと球面S2との形状の相違の程度を評価することが可能となる。ここで、式61に示したベクトルδSD2を最小化するため、式62に示すような、目的関数歪評価関数gSD2eを導入する。目的関数gSD2eを最小化する解e=e02は、以下の式63式で表現される連立方程式に対して、ガウス消去法や、コレスキー分解法などを適用することで求められる。なお、式63は、式62に対して、最小二乗法を適用することによって算出される。解e02に基づき、式38により歪評価行列Eを復元し、式43に示した歪評価値gDEの値、即ち、評価行列Eのノルムを求める。そして、歪評価値gDEの値が、閾値δ0以下であるか否かを判定する。判定結果が肯定の場合、CPU100は、処理を後続のオフセット更新処理に進める。一方、判定結果が否定である場合には、CPU100は、一連の処理を終了する。[9.オフセット更新処理]オフセット更新処理は、歪評価値gDEの値が、閾値δ0以下である場合に実行される処理である。オフセット更新処理では、中心点算出処理において算出された中心点x0の座標が示す値を新オフセットc1として採用したうえで、オフセットcを、旧オフセットc0から、新オフセットc1へと更新する。本実施形態では、中心点x0を新オフセットc1として採用しているが、球面S2の中心点x02を新オフセットc1として採用しても良い。歪評価値gDEの値が閾値δ0以下である場合には、球面Sの中心点x0と、球面S2の中心点x02とは、いずれも、地磁気Bgを表す球面SGの中心点x0Gとほぼ一致するからである。なお、球面S2の中心点x02は、目的関数gSD2eを最小化する解e=e02のうち、式51のeXに相当する部分を、式52に代入することにより求められる。[10.結論]以上に示したように、本実施形態は、複数の磁気データq1~qNの分布の3次元的な広がりの程度を表す共分散行列Aの最小固有値λ3と広がりの最も小さな方向を表す共分散行列Aの最小固有値λ3に対応する固有ベクトルu3とを算出する磁気データ分布指標算出処理、複数の磁気データq1~qNの分布が3次元的な場合及び2次元的な場合の双方において新オフセットc1の候補である中心点x0を算出する中心点算出処理、複数の磁気データq1~qNの分布が3次元的な場合及び2次元的な場合の双方において立体SDの形状の歪の程度を評価して中心点x0をオフセットとして採用することが妥当であるかを判断する歪判定処理、及び歪判定処理における判定結果が肯定である場合に中心点x0を新オフセットc1に採用してオフセットcを旧オフセットc0から新オフセットc1へと更新するオフセット更新処理を有する。これらの4つの処理により、複数の磁気データq1~qNの分布が3次元的な場合及び2次元的な場合の双方において、外部磁界Bxの影響が大きい環境下で測定された複数の磁気データq1~qNに基づいて算出された中心点x0をオフセットcとして採用することを防止することで、内部磁界Biに近い値を示す適切な中心点x0のみをオフセットcとして採用することを可能とした。すなわち、本実施形態は、複数の磁気データq1~qNの分布が3次元的な場合及び2次元的な場合の双方において、正確なオフセットによる、正確な地磁気Bgの値の算出が可能になった。これにより、3次元磁気センサ60が、自動車等のように左右方向のみに姿勢を変化させ、上下方向の姿勢変化が小さい機器1に搭載されている場合等にも、適切な地磁気Bgを算出することが可能となった。歪判定処理では、複数の磁気データq1~qNの座標と立体SDの表面との誤差を表すベクトルである誤差δSD、旧オフセットc0から球面S2の中心点x02を示すベクトルである更新ベクトルk2と複数の磁気データq1~qNの分布の広がりが最も小さい方向を示す固有ベクトルu3との内積を表す値である誤差δU2、及び複数の磁気データq1~qNの分布が2次元的な広がりを有する場合に広がりを有する方向のみについて立体SDと球面S2との形状の相違の程度を評価することを表す誤差δEを要素とするベクトルδSD2を定義した。そして、ベクトルδSD2のノルムを最小化するような歪評価行列E及び中心点x02を算出し、歪評価行列Eの最大固有値λE1の大きさ歪評価行列Eのノルムを表す歪評価値gDEに基づいて、立体SDと球面S2との形状の相違の程度を評価した。複数の磁気データq1~qNが2次元的に分布する場合には、式53に示される、9行9列の行列X2TX2の階数が、9よりも小さくなる場合がある。この場合には、式53に基づいて、解e0を算出することはできず、評価行列Eのノルムを表す歪評価値gDEを算出することができないため、立体SDと球面S2との形状の相違の程度を評価することができなかった。それに対して、本実施形態では、誤差δU2、誤差δEを導入することで、複数の磁気データq1~qNが2次元的に分布する場合も考慮された式63によって解e02を算出し、解e02に基づいて歪評価値gDEを算出した。これにより、本実施形態に係る歪判定処理は、複数の磁気データq1~qNが2次元的に分布する場合であっても、立体SDと球面S2との形状の相違の程度を評価することが可能となった。そして、外部磁界Bxの影響を受けた複数の磁気データq1~qNに基づいて算出された不適切な中心点x0をオフセットcとして採用することを防止し、内部磁界Biに近い値を示す適切な中心点x0のみをオフセットcとして採用することを可能とした。また、歪評価行列Eは、複数の磁気データq1~qNの分布が2次元的な広がりを有する場合に固有ベクトルu3方向について立体SDと球面S2との形状の相違の程度を評価しないことを表す誤差δEを要素とするベクトルδSD2のノルムを最小化することで算出される。歪評価行列Eを算出するにあたって、誤差δEを最小化することにより、複数の磁気データq1~qNの分布が2次元的な場合に、立体SDと球面S2との形状の相違の程度について評価を行うことができない固有ベクトルu3方向については、立体SDと球面S2との形状の相違の程度について評価を行わず、固有ベクトルu1及びu2方向についてのみ、立体SDと球面S2との形状の相違の程度を評価する歪評価行列Eを算出することが可能となった。なお、誤差δEは、固有ベクトルu3方向について立体SDと球面S2との形状の相違の程度を評価しないことを表す第2の制約条件式54参照と、第2の制約条件を複数の磁気データq1~qNの3次元的な広がりの程度固有ベクトルu3方向の広がりの程度を表す固有値λ3の大きさに応じて緩和する係数βとを、掛け算したものである。すなわち、本実施形態に係る歪判定処理は、このような係数βを有する誤差δEを考慮して歪評価行列Eを算出するものであるため、複数の磁気データq1~qNが2次元的に分布している場合のみならず、3次元的な広がりを有して分布している場合についても、3次元的な広がりの程度に応じて、立体SDと球面S2との形状の相違の程度を評価することが可能となった。また、中心点x0は、誤差δUを要素とするベクトルであるベクトルδSUのノルムを最小化することで算出される。中心点x0を算出するにあたって、誤差δUを最小化することにより、複数の磁気データq1~qNの分布が2次元的な場合に、球面Sπの中心点xπが有する情報のうち、固有ベクトルu3方向の不正確な情報を排除して、新オフセットc1の候補である球面Sの中心点x0を算出することが可能となる。これにより、本実施形態に係る中心点算出処理は、複数の磁気データq1~qNが2次元的に分布する場合にも、3次元的な広がりの程度に応じて、正確なオフセットcの候補である球面Sの中心点x0を算出することが可能となった。また、本実施形態に係る磁気データ分布指標算出処理は、複数の磁気データq1~qNのセンサ座標系ΣSにおける分布の3次元的な広がりの程度を表す、共分散行列Aの最小固有値λ3と、広がりが最も狭い方向を表す、固有ベクトルu3とを算出する。これにより、中心点算出処理において、複数の磁気データq1~qNの分布の3次元的な広がりの程度に応じた、適切な中心点x0を算出することを可能とすると共に、歪判定処理において、複数の磁気データq1~qNの分布の3次元的な広がりの程度に応じて、立体SDと球面S2との形状の相違の程度を判定することが可能となった。なお、本実施形態に係る歪判定処理は、式41に示した誤差δSDを要素とするベクトルδSD2のノルムを最小化することで算出される歪評価行列Eに基づいて、立体SDと球面S2との形状の相違の程度を評価する。誤差δSDは、複数の磁気データq1~qNが、ある球面S2から、単にホワイトノイズに相当する誤差を有して分布する場合と、球面S2とは異なる歪んだ形状の立体SDの表面近傍に分布する場合とを、区別して表現する。例えば、式17に示した誤差δSは、N個の磁気データq1~qNが球面Sから誤差を有して分布する場合に、その誤差が、ホワイトノイズによるものであるのか、それとも、不均一な外部磁界Bxによるものであるのかを判断できなかった。つまり、N個の磁気データq1~qNが球面S2から大きく歪んだ形状を有する立体SDの表面近傍に分布している場合であっても、単純に球面からのホワイトノイズとしての誤差として評価するにとどまっていた。これに対して、本実施形態に係る歪判定処理は、誤差δSDを用いるため、ホワイトノイズによる誤差と、不均一な外部磁界Bxの影響による誤差とを区別することが可能となった。N個の磁気データq1~qNが示すN個の座標と球面S2とのN個の誤差が、全体として単なるホワイトノイズである場合、中心点算出処理において算出された中心点x0と、内部磁界Biを表すベクトル地磁気Bgを表す球面SGの中心点x0Gとは、ほぼ同一の座標であると看做すことができる。本実施形態における歪判定処理は、誤差が単なるホワイトノイズである場合と、不均一な外部磁界Bxによる歪であるかを判別することが可能であるため、単なるホワイトノイズのみが存在する場合に算出された中心点x0を、オフセットcとして採用することが可能である。一方、本実施形態における歪判定処理は、不均一な外部磁界Bxによる立体SDの歪が存在する場合に算出された中心点x0を、オフセットcとして採用することを防止することができる。これにより、本実施形態は、不正確な値をオフセットとして採用することを未然に防止することが可能になるという利点を有する。なお、歪評価行列Eは、N個の磁気データq1~qNが示すN個の座標と球面S2とのN個の誤差の大きさを表す行列であり、中心点x0から見て、磁気データqiの示す座標と球面S2との誤差が大きな方向を向く固有ベクトルuE1と、誤差の大きさを表す最大固有値λE1とを有する。本実施形態の歪判定処理は、歪評価行列Eの最大固有値λE1の大きさ歪評価行列Eのノルムを表す歪評価値gDEに基づいて、立体SDと球面S2との形状の相違の程度を評価することを可能とした。これにより、外部磁界Bxの影響が大きく、立体SDと球面S2との形状の相違の程度が大きい場合に、中心点算出処理において算出される球面Sの中心点x0をオフセットとして採用することを防止することが可能になった。また、本実施形態に係る中心点算出処理は、複数の磁気データq1~qNが、中心点x0を中心とした球面Sの近傍に分布すると仮定することで、複数の磁気データq1~qNより中心点x0を算出し、これをオフセットの候補とした。このような中心点算出処理は、オフセットの候補を簡単な演算により算出できるという利点を有し、処理速度の向上という利点を有するものである。また、本実施形態に係る中心点算出処理は、中心点x0の算出にあたり、複数の磁気データq1~qNの重心qcを原点とする重心座標系ΣCにおいて、複数の磁気データq1~qNの値を表現することで計算を行った。これにより、3次元磁気センサ60の出力値である、センサ座標系ΣSにおける複数の磁気データq1~qNの値をそのまま用いて計算する場合に比べて、データの値を小さくすることができ、メモリサイズの小さい変数型で表現することが可能となった。すなわち、本実施形態の磁気データ処理プログラム70は、メモリ使用量の低減化を可能とするという利点を有する。<B.変形例>本発明は上述した実施形態に限定されるものではなく、例えば、以下の変形が可能である。また、以下に示す変形例のうちの2以上の変形例を組み合わせることもできる。1変形例1上述した実施形態では、係数αを、式26に示すような、固有値λ3の関数として定められているが、係数αは2値で与えられても良い。例えば、固有値λ3が閾値λ0以上となる場合には、係数αを0または0に限りなく近い正の小さな値第2の値となるように設定し、固有値λ3が閾値λ0未満となる場合には、係数αを正の大きな値第1の値となるように設定しても良い。この場合、固有値λ3が閾値λ0未満の値となり、中心点xπの座標が有する情報が固有ベクトルu3方向に正確性を有さないと判断される場合に、更新ベクトルk及び中心点x0と更新ベクトルk2及び中心点x02とを、中心点xπの座標が有する情報のうち正確性の高い固有ベクトルu1及びu2方向の情報のみを利用して生成することが可能となる。同様に、上述した実施形態では、係数βを、式56に示すような、固有値λ3の関数として定められているが、係数βは2値で与えられても良い。例えば、固有値λ3が閾値λ0以上となる場合には、係数βを0または0に限りなく近い正の小さな値第4の値となるように設定し、固有値λ3が閾値λ0未満となる場合には、係数βを正の大きな値第3の値となるように設定しても良い。この場合、固有値λ3が閾値λ0未満の値となり、複数の磁気データq1~qNの分布が2次元的な場合に、立体SDと球面S2との形状の相違の程度をより正確に評価することが可能となる。2変形例2上述した実施形態では、係数αを、式26に示すような、固有値λ3の関数として定められているが、係数αを固有値λ1及び固有値λ3の大きさの比率を変数とする関数としても良い。例えば、固有値λ3が固有値λ1に比べて小さくなるに従って、係数αの値が大きくなるような関数としても良い。同様に、上述した実施形態では、係数βを、式56に示すような、固有値λ3の関数として定められているが、係数βを固有値λ1及び固有値λ3の大きさの比率を変数とする関数としても良い。例えば、固有値λ3が固有値λ1に比べて小さくなるに従って、係数βの値が大きくなるような関数としても良い。この場合も、複数の磁気データq1~qNの分布が2次元的な場合に、立体SDと球面との形状の相違の程度をより正確に評価することが可能であり、また、複数の磁気データq1~qNの分布が2次元的な場合に、中心点xπが有する情報のうち固有ベクトルu3方向の不正確な情報を排除して新オフセットc1の候補である中心点x0を算出することが可能である。3変形例3上述した実施形態及び変形例では、式22に示す制約条件を係数αにより緩和した誤差δU2と、式54に示す制約条件を係数βにより緩和した誤差δEとを要素として含む、ベクトルδSD2の大きさを表す目的関数gSD2eを最小化する解e02に基づいて、歪評価行列Eを算出しているが、本発明はこのような形態に限定されるものでは無く、共分散行列Aの最小固有値λ3が閾値λ0未満の値となり複数の磁気データq1~qNの分布が2次元的な広がりを有する3次元的な広がりを有さないと看做せる場合に、式22及び式54に示す制約条件を厳密に満たすという条件の下で、歪評価行列Eを算出しても良い。例えば、複数の磁気データq1~qNの分布が2次元的な場合、式22及び式54に示す制約条件を厳密に満たすという条件の下で、式50に示した目的関数gSDeを最小化するような解e=e0Aを算出し、この解を用いて歪評価行列Eを算出しても良い。式22及び式54に示す制約条件を厳密に満たすという条件の下で算出された、目的関数gSDeを最小化するような解e0Aは、ラグランジュ未定乗数法を用いて以下の式64により算出される。なお、前述の通り、式22に示す制約条件は、厳密には、更新ベクトルkと固有ベクトルu3とが直交する条件を表しているが、本変形例では、8.2節と同様に、式22に示す制約条件を、更新ベクトルk2と固有ベクトルu3とが直交することを表す制約条件として用いる。ここで、行列M2は、式65に示す4行9列の行列であり、ベクトルρは、4つのラグランジュ未定乗数を要素とする、4次元のベクトルである。また、ベクトルbρは、式66に示す4次元のベクトルである。式64により算出される解e0Aは、式22及び式54に示す制約条件を厳密に満たす解であるため、複数の磁気データq1~qNの分布が2次元的な場合に、立体SDと球面との形状の相違の程度をより正確に評価することが可能となる。4変形例4上述した実施形態及び変形例では、式22に示す制約条件を係数αにより緩和した誤差δUを要素として含むベクトルであるベクトルδSUを最小化することで式35により中心点x0を算出しているが、本発明はこのような形態に限定されるものでは無く、共分散行列Aの最小固有値λ3が閾値λ0未満の値となり、複数の磁気データq1~qNの分布が2次元的な広がりを有する3次元的な広がりを有さないと看做せる場合に、式22に示す制約条件を厳密に満たすような中心点xAを算出しても良い。式22に示す制約条件を厳密に満たすような中心点xAは、ラグランジュ未定乗数法を用いて以下の式67により算出される。ここで、行列X4は、式68に示す4行4列の行列であり、ベクトルx4は、ラグランジュ未定乗数ρ4と、式18に示す3次元の変数を表すベクトルxとを要素とする、4次元のベクトルである。また、ベクトルb4は、式70に示す4次元のベクトルである。式67により算出される中心点xAは、式22に示す制約条件を厳密に満たす座標であるため、複数の磁気データq1~qNの分布が2次元的な場合に、中心点xπが有する情報のうち固有ベクトルu3方向の不正確な情報を厳密に排除して新オフセットc1の候補である中心点xAを算出することが可能となる。5変形例5上述した実施形態及び変形例では、中心点算出処理において、式22に示す制約条件を係数αにより緩和した誤差δUを要素として含むベクトルであるベクトルδSUを最小化することで式35により中心点x0を算出しているが、本発明はこのような形態に限定されるものでは無く、式22に示す制約条件を考慮せずに中心点x0を算出するものであっても良い。例えば、2節2.磁気データが球面近傍に分布する場合の球面の中心点の算出に示した方法により、式20を用いて中心点x0を求めても良い。また、磁気データ分布指標算出処理において算出される共分散行列Aの最小固有値λ3が閾値λ0以上の値となり、複数の磁気データq1~qNの分布が3次元的な広がりを有すると看做せる場合にのみ、式20を用いて中心点x0を求めても良い。中心点x0は、複数の磁気データq1~qNが、ある球面の近傍に分布することを仮定した場合の当該球面の中心点である。従って、中心点x0は、複数の磁気データq1~qNを近傍に有する球面の中心点を算出可能な任意の方法によって算出することができる。そして、中心点x0の算出に複数の磁気データq1~qNを用いる場合、複数の磁気データq1~qNを表面近傍に有する立体SDの形状と球面の形状との相違の程度を判定することにより、中心点x0が不均一な外部磁界Bxの影響下で算出された誤差を含む値であるか否か、すなわち、中心点x0をオフセットcとして採用することが適当であるか否かを判断することが可能になる。6変形例6上述した実施形態及び変形例では、誤差δSD、誤差δU2、及び誤差δEを要素とするベクトルδSD2のノルムを最小化するような歪評価行列Eを算出することにより、複数の磁気データq1~qNが2次元的に分布している場合及び3次元的に分布している場合の双方の場合を想定し、複数の磁気データq1~qNの3次元的な広がりの程度に応じて、立体SDと球面との形状の相違の程度を評価するものであるが、本発明はこのような形態に限定されるものでは無く、複数の磁気データq1~qNが3次元的な広がりを有して分布していると仮定して、立体SDと球面との形状の相違の程度を評価しても良い。例えば、8.1節8.1.磁気データが3次元的に分布する場合の歪判定処理において述べた、複数の磁気データq1~qNが3次元的な広がりを有して分布する場合に限定した歪判定処理により、立体SDと球面との形状の相違の程度を評価しても良い。具体的には、式50で示した目的関数gSDeを最小化する解e=e0を式53により算出したうえで、解e0を用いて算出される歪評価行列Eのノルムを求めることで、立体SDと球面との形状の相違の程度を評価することができる。また、磁気データ分布指標算出処理において算出される共分散行列Aの最小固有値λ3が閾値λ0以上の値となり、複数の磁気データq1~qNの分布が3次元的な広がりを有すると看做せる場合にのみ、3次元的な広がりを有して分布する場合に限定した歪判定処理を用いて、立体SDと球面との形状の相違の程度を評価しても良い。この場合、磁気データ分布指標算出処理において算出される共分散行列Aの最小固有値λ3が閾値λ0未満の値であれば、初期化処理及び磁気データ取得処理を再度実施しても良い。すなわち、この変形例から、以下の発明が把握される。磁界を発生させる部品を備えた機器に組み込まれ、互いに直交する3方向の磁気成分をそれぞれ検出する3次元磁気センサと、前記3次元磁気センサから、3軸の座標系において表現される3次元のベクトルデータとして順次出力される磁気データを蓄積する蓄積手段と、前記蓄積手段に蓄積された前記複数の磁気データで示される複数の座標の分布の3次元的な広がりの程度を示す分散評価値を算出する、分布指標算出手段と、前記複数の磁気データの各々が示す3軸の座標が、中心点算出用球面の近傍に確率的に分布すると仮定して、前記複数の磁気データに基づいて前記中心点算出用球面の中心点を算出する、中心点算出手段と、前記複数の磁気データの各々が示す3軸の座標が、歪判定用球面を歪ませた形状の立体の表面近傍に確率的に分布すると仮定して、前記複数の磁気データに基づいて、前記立体と、前記歪判定用球面との形状の相違の程度を示す歪評価値を算出し、前記歪評価値が歪許容値以下であるか否かを判定する、歪判定手段と、前記歪判定手段の判定結果が肯定である場合に、前記磁気データのうち前記部品の発生する磁界の成分を表す3軸の座標であるオフセットを、旧オフセットから前記中心点に更新する、オフセット更新手段と、を備えることを特徴とする、地磁気測定装置。7変形例7上述した実施形態及び変形例において、オフセットcを算出し更新する処理は、初期化処理、磁気データ取得処理、磁気データ分布指標算出処理、中心点算出処理、歪判定処理、及びオフセット更新処理を含むが、本発明はこのような形態に限定されるものでは無く、磁気データ分布指標算出処理を備えずに構成されても良い。この場合、中心点算出処理は、中心点x0を、式20を用いて算出するものであっても良い。また、この場合、立体SDと球面との形状の相違の程度を評価は、8.1節に示した複数の磁気データq1~qNが3次元的な広がりを有して分布する場合に限定した歪判定処理により行い、歪評価行列Eを式53により求められる解e0を用いて算出しても良い。8変形例8上述した実施形態及び変形例において、オフセットcを算出し更新する処理は、初期化処理、磁気データ取得処理、中心点算出処理、歪判定処理、及びオフセット更新処理を含むが、本発明はこのような形態に限定されるものでは無く、中心点算出処理を備えずに構成されても良い。9変形例9上述した実施形態及び変形例において、歪を表す成分kEは、式37に表したように、式39に示した関数kepに対して、球面Sの中心点x0を原点とする座標系において磁気データqiの座標を示すベクトルqi-x0を代入することにより与えられる値を各要素とするベクトルであるが、本発明はこのような形態に限定されるものでは無く、歪を表す成分kEは、球面Sの中心点x0とは異なる点を原点とする座標系において磁気データqiの座標を示すベクトルを式39に代入して得られる値を各要素としても良い。例えば、歪を表す成分kEは、以下の式71のように、複数の磁気データq1~qNの重心qcを原点とする座標系において磁気データqiの座標を示すベクトルを用いて表されるものであって良い。また、前述の中心点算出処理よりも計算量の少ない方法により、複数の磁気データq1~qNを近傍に有する球面の中心点の候補となる座標を求め、当該座標を原点とする座標系において磁気データqiの座標を示すベクトルを用いて歪を表す成分kEを表しても良い。中心点算出手段により算出される球面Sの中心点x0とは異なる座標を原点とする座標系において磁気データqiの座標を示すベクトルを用いて歪を表す成分kEを表現する場合であっても、歪評価行列Eの各成分は、歪評価行列Eの最大固有値λE1に対応する固有ベクトルuE1の方向と、球面S2からの誤差が大きな磁気データqiが多く存在する位置を当該座標から見た方向とが近づくように定められる。また、歪評価行列Eの最大固有値λE1は、当該座標から見て、固有ベクトルuE1の方向に存在する複数の磁気データと球面S2との誤差の大きさの程度を表す値となる。従って、この場合も、歪評価行列Eによって、立体SDと球面S2との形状の相違の程度を評価することが可能となる。なお、球面Sの中心点x0以外の座標を原点とする座標系において磁気データqiの座標を示すベクトルを用いて歪を表す成分kEを表現する場合、中心点算出処理は、歪判定処理による判定結果が肯定である場合歪が小さい場合のみ実施しても良い。この場合、中心点算出処理を行うことなく歪判定処理を行うことが可能となる。すなわち、中心点算出処理よりも計算量の小さな方法により、複数の磁気データq1~qNを近傍に有する球面の中心点の候補となる座標を求め、当該座標を用いて歪判定処理を行うため、歪判定処理による判定結果が否定である場合の計算負荷を低減することが可能となる。また、球面Sの中心点x0以外の座標を原点とする座標系において磁気データqiを示すベクトルを用いて歪を表す成分kEを表現する場合、中心点算出処理を実施せず、歪判定処理において算出される球面S2の中心点x02を新オフセットc1として採用しても良い。この場合、中心点算出処理を実施することなく、オフセットcを算出し更新する処理を行うことが可能になるため、オフセットcを算出し更新する処理の簡素化と、計算負荷の低減による処理の高速化が可能となる。10変形例10上述した実施形態及び変形例において、歪判定処理は、球面Sの中心点x0を原点とする座標系において磁気データqiの座標を示すベクトルqi-x0を、関数kepに対して代入することにより与えられる値を各要素とするベクトルである歪を表す成分kEを用いて算出された歪評価行列Eにより、立体SDと球面S2との形状の相違の程度を評価するものであるが、本発明はこのような形態に限定されるものでは無く、球面Sの中心点x0を原点とする座標系において磁気データqiの座標を示すベクトルを関数kepに対して代入して得られる値を要素とする歪を表す成分kEを用いて算出される歪評価行列Eと、球面Sの中心点x0とは異なる座標を原点とする座標系において磁気データqiの座標を示すベクトルを関数kepに対して代入して得られる値を要素とする歪を表す成分kEを用いて算出される歪評価行列Eとの双方を利用して、立体SDと球面S2との形状の相違の程度を評価しても良い。例えば、歪判定処理は、まず、重心qcを原点とする座標系において磁気データqiの座標を示すベクトルにより表される歪を表す成分kEを用いて算出される歪評価行列Eにより、立体SDと球面S2との形状の相違の程度を評価し第1の評価、次に、第1の評価において、立体SDと球面S2との形状の相違の程度が小さいと評価された場合に、中心点x0を原点とする座標系において磁気データqiの座標を示すベクトルにより表される歪を表す成分kEを用いて算出される歪評価行列Eにより、立体SDと球面S2との形状の相違の程度について評価第2の評価するものであっても良い。そして、第1の評価及び第2の評価の双方において、立体SDと球面S2との形状の相違の程度が小さいと評価された場合に、歪判定処理による判定結果は肯定形状の相違の程度は小さいと判定されることとなり、中心点x0を新オフセットc1として採用するものであっても良い。また、歪判定処理は、まず、第2の評価を行い、次に、第2の評価において立体SDと球面S2との形状の相違の程度が小さいと評価された場合に、第1の評価を実施するものであっても良い。このように、異なる2つの座標系において表現された2つの歪を表す成分kEを用いて、立体SDと球面S2との形状の相違の評価を2回行うことにより、1つの座標系において表現された1つの歪を表す成分kEを用いて1回のみ評価を行う場合に比べて、立体SDと球面S2との形状の相違をより正確に評価することが可能となる。11変形例11また、上述した実施形態では、機器1は、CPU10、RAM20、ROM30、通信部40、表示部50、及び3次元磁気センサ60を備えるが、本発明はこのような形態に限定されるものでは無く、機器1は、3次元磁気センサ60以外のセンサを備えるものであっても良い。例えば、機器1が、3次元磁気センサ60の他に、加速度センサや角速度センサ等の3次元磁気センサ60以外のセンサを備える場合、これら複数のセンサからの出力を用いることにより、内部磁界Biの方向及び大きさを表すベクトルSBi算出し、これを中心点x0としても良い。歪判定処理は、どのような方法により中心点x0を算出した場合であっても、当該中心点x0が不均一な外部磁界Bxの影響を受けないで算出された、内部磁界Biの方向及び大きさを適切に表すものであるか否かを判定することが可能である。"