"軟磁性材料を有する部品を備えた機器に組み込まれ、互いに直交する第1軸、第2軸、及び第3軸の3方向の磁気成分をそれぞれ検出する3次元磁気センサと、前記3次元磁気センサから、3軸の座標系において表現される3次元のベクトルデータとして順次出力される磁気データを蓄積する蓄積部と、複数の前記磁気データの各々が示す3軸の座標を近傍に有し、互いに形状の異なる第1楕円面、第2楕円面、及び第3楕円面のうち、少なくとも2つの楕円面の各々の中心点を表す3軸の座標を算出する初期楕円面生成部と、前記少なくとも2つの楕円面の各々の中心点を表す3軸の座標の相互間の距離が、第1閾値以下であるか否かを判定する、初期楕円面中心点判定部と、前記初期楕円面中心点判定部の判定結果が肯定である場合、前記少なくとも2つの楕円面のうち、少なくとも1つの楕円面の形状を表す係数行列に基づいて、初期楕円面補正行列を算出すると共に、前記少なくとも1つの楕円面の中心点を表す3軸の座標に基づいて、初期中心点の3軸の座標を算出する、初期補正値生成部と、を備えることを特徴とする地磁気測定装置。前記初期楕円面生成部は、複数の前記磁気データの各々が示す3軸の座標が、楕円面の近傍に確率的に分布すると仮定して、前記楕円面を、前記第1軸成分の二乗を表す項、前記第2軸成分の二乗を表す項、及び前記第3軸成分の二乗を表す項を含む楕円面方程式により表現したとき、複数の前記磁気データの各々が表す3軸の座標を前記楕円面方程式のうち前記第1軸成分の二乗を表す項以外の項に代入して得られる値と、複数の前記磁気データの各々が表す3軸の座標のうち前記第1軸成分の二乗値との誤差を最小化する前記第1楕円面の中心点を表す3軸の座標を算出する、第1楕円面生成部、複数の前記磁気データの各々が表す3軸の座標を前記楕円面方程式のうち前記第2軸成分の二乗を表す項以外の項に代入して得られる値と、複数の前記磁気データの各々が表す3軸の座標のうち前記第2軸成分の二乗値との誤差を最小化する前記第2楕円面の中心点を表す3軸の座標を算出する、第2楕円面生成部、及び、複数の前記磁気データの各々が表す3軸の座標を前記楕円面方程式のうち前記第3軸成分の二乗を表す項以外の項に代入して得られる値と、複数の前記磁気データの各々が表す3軸の座標のうち前記第3軸成分の二乗値との誤差を最小化する前記第3楕円面の中心点を表す3軸の座標を算出する、第3楕円面生成部のうち、少なくとも2つを備えることを特徴とする、請求項1に記載の地磁気測定装置。3次元の変数ベクトルの示す3軸の座標を起点として前記磁気データの示す3軸の座標を表す3次元のベクトルを第1変数ベクトルとし、前記第1変数ベクトルを、変数行列により変換して得られる3次元のベクトルを第2変数ベクトルとし、前記第2変数ベクトルの示す3軸の座標を表すデータを変換後データとし、複数の前記変換後データの各々が示す3軸の座標と、前記変数ベクトルの示す3軸の座標を中心とする球面との誤差を表し、前記変数行列の各成分と、前記変数ベクトルの各要素とを変数とする関数を、楕円面最適化関数としたとき、前記楕円面最適化関数の変数の初期値として、前記初期楕円面補正行列の各成分と、前記初期中心点の示す3軸の座標とを適用し、その後、前記楕円面最適化関数の変数の値を逐次更新することにより、前記楕円面最適化関数を最小化する解として、最適楕円面補正行列と、最適中心点の示す3軸の座標とを算出する、楕円面最適補正値生成部と、前記3次元磁気センサから出力される磁気データの示す3軸の座標を、前記最適中心点の示す3軸の座標から表した3次元のベクトルを、前記最適楕円面補正行列により変換することにより、地磁気の向きを算出する、地磁気算出部と、を更に備えることを特徴とする、請求項1または2に記載の地磁気測定装置。前記3次元磁気センサから出力される磁気データの示す3軸の座標を、前記初期中心点の示す3軸の座標から表した3次元のベクトルを、前記初期楕円面補正行列により変換することにより、地磁気の向きを算出する、地磁気算出部を更に備えることを特徴とする、請求項1または2に記載の地磁気測定装置。しかし、非特許文献1及び非特許文献2に開示された方法で算出される座標変換行列の初期値、つまり、座標変換行列の候補となる行列は、非線形関数の大局的最適解となる座標変換行列とは大きく異なる値である場合がある。非線形最適化演算に用いられる初期値が、非線形関数の大局的最適解と大きく異なる値である場合、非線形最適化演算により求められる最適解が、大局的最適解とは異なる局所最適解となる可能性が高い。従って、非特許文献1及び非特許文献2に開示された方法で算出される座標変換行列を用いて、3次元磁気センサから出力されるベクトルデータを座標変換しても、正確な地磁気の向きを求めることはできない可能性が高いという問題点があった。本発明は、上述した事情を鑑み、非線形最適化問題の対象となる非線形関数の大局的最適解に近い正確な初期値を算出し、当該初期値に基づいて座標変換行列を算出することにより、正確な地磁気の方向を算出することを解決課題とする。以下、本発明について説明する。なお、本発明の理解を容易にするために本実施形態、変形例、及び添付図面の参照符号を括弧書きにて付記するが、それにより本発明が本実施形態に限定されるものではない。上述した課題を解決するため、本発明に係る地磁気測定装置は、軟磁性材料を有する部品を備えた機器に組み込まれ、互いに直交する第1軸、第2軸、及び第3軸の3方向の磁気成分をそれぞれ検出する3次元磁気センサと、前記3次元磁気センサから、3軸の座標系において表現される3次元のベクトルデータとして順次出力される磁気データを蓄積する蓄積部と、複数の前記磁気データの各々が示す3軸の座標を近傍に有し、互いに形状の異なる第1楕円面、第2楕円面、及び第3楕円面のうち、少なくとも2つの楕円面の各々の中心点を表す3軸の座標を算出する初期楕円面生成部と、前記少なくとも2つの楕円面の各々の中心点を表す3軸の座標の相互間の距離が、第1閾値以下であるか否かを判定する、初期楕円面中心点判定部と、前記初期楕円面中心点判定部の判定結果が肯定である場合、前記少なくとも2つの楕円面のうち、少なくとも1つの楕円面の形状を表す係数行列に基づいて、初期楕円面補正行列を算出すると共に、前記少なくとも1つの楕円面の中心点を表す3軸の座標に基づいて、初期中心点の3軸の座標を算出する、初期補正値生成部と、を備えることを特徴とする。ソフトアイアン効果が生じている場合、複数の磁気データの示す座標は、楕円面近傍に分布するため、地磁気の向きを算出するためには、楕円面近傍に分布する座標を、当該楕円面と同一の中心点を有する球面近傍に分布する座標へと変換する座標変換つまり、楕円補正を行う必要がある。このような座標変換を行うための行列を算出するためには、複数の磁気データを近傍に有し、複数の磁気データの分布形状を正確に表現する楕円面の形状を特定することが必要となる。しかし、複数の磁気データの示す座標の分布の形状から、楕円面の形状を特定することが困難な場合、例えば、複数の磁気データの示す座標の分布が広がりを有さず特定の座標の近傍に分布する場合であっても、複数の磁気データの示す座標を近傍に有するような楕円面を強引に算出することができる。このような楕円面は、複数の磁気データの分布形状を正確に表現したものではない。複数の磁気データの示す座標の分布形状を正確に表さない不適切な楕円面に基づいて生成される座標変換行列により楕円補正を行う場合、正確な地磁気の向きを算出することは困難である。従って、複数の磁気データの示す座標の分布の形状から、楕円面の形状を特定することが困難な場合には、楕円補正を行うための座標変換行列が算出されることを防止する必要がある。この発明によれば、初期楕円面生成部は、複数の磁気データの示す座標を近傍に有する楕円面として、第1楕円面、第2楕円面、及び、第3楕円面の3つの異なる楕円面のうち、少なくとも2つの楕円面の各々の中心点の座標を算出する。そして、初期楕円面中心点判定部は、2つの楕円面各々の中心点の相互間の距離が、第1閾値以下であることを判定する。初期楕円面中心点判定部の判定結果が肯定である場合、初期楕円面生成部が生成する少なくとも2つの楕円面は、いずれも、複数の磁気データの示す座標を近傍に有し、初期楕円面生成部が生成する少なくとも2つの楕円面各々の中心点の示す座標は、同一と看做すことができる。従って、初期楕円面中心点判定部の判定結果が肯定である場合、初期楕円面生成部が生成する少なくとも2つの楕円面は、同一の形状を有すると看做すことができる。複数の磁気データの示す座標が、楕円面の形状を特定することができるような形状に分布している場合、初期楕円面生成部において生成される少なくとも2つの異なる楕円面は、複数の磁気データの分布により特定される楕円面の形状と、同一と看做すことのできる形状を有する楕円面として、各々算出される。逆に、複数の磁気データの示す座標の分布形状から、楕円面の形状を特定することが困難な場合、初期楕円面生成部が生成する少なくとも2つの楕円面の形状の各々は、複数の磁気データの示す座標を近傍に有するという条件のみに基づいて定められる。この場合、少なくとも2つの楕円面の形状は、同一の形状と看做すことができない可能性が高く、少なくとも2つの楕円面の中心点の各々が示す座標も、同一の座標とは看做すことができない可能性が高い。本発明に係る地磁気測定装置は、第1楕円面、第2楕円面、及び第3楕円面のうち、少なくとも2つの楕円面が同一の形状と看做すことができることを判定したうえで、初期楕円面補正行列を生成する。従って、複数の磁気データの示す座標の分布の形状から、楕円面の形状を特定することが困難な場合に、不適切な初期楕円面補正行列が生成されることを防止することができる。また、上述した地磁気測定装置において、前記初期楕円面生成部は、複数の前記磁気データの各々が示す3軸の座標が、楕円面の近傍に確率的に分布すると仮定して、前記楕円面を、前記第1軸成分の二乗を表す項、前記第2軸成分の二乗を表す項、及び前記第3軸成分の二乗を表す項を含む楕円面方程式により表現したとき、複数の前記磁気データの各々が表す3軸の座標を前記楕円面方程式のうち前記第1軸成分の二乗を表す項以外の項に代入して得られる値と、複数の前記磁気データの各々が表す3軸の座標のうち前記第1軸成分の二乗値との誤差を最小化する前記第1楕円面の中心点を表す3軸の座標を算出する、第1楕円面生成部、複数の前記磁気データの各々が表す3軸の座標を前記楕円面方程式のうち前記第2軸成分の二乗を表す項以外の項に代入して得られる値と、複数の前記磁気データの各々が表す3軸の座標のうち前記第2軸成分の二乗値との誤差を最小化する前記第2楕円面の中心点を表す3軸の座標を算出する、第2楕円面生成部、及び、複数の前記磁気データの各々が表す3軸の座標を前記楕円面方程式のうち前記第3軸成分の二乗を表す項以外の項に代入して得られる値と、複数の前記磁気データの各々が表す3軸の座標のうち前記第3軸成分の二乗値との誤差を最小化する前記第3楕円面の中心点を表す3軸の座標を算出する、第3楕円面生成部のうち、少なくとも2つを備えることを特徴とする、ことが好ましい。複数の磁気データの示す座標が、楕円面の形状を特定することができるような形状に分布している場合、複数の磁気データの示す座標と楕円面との誤差をどのような形式により表現しても、複数の磁気データの示す座標と楕円面との誤差を最小化することにより生成される楕円面は、複数の磁気データの分布により特定される楕円面の形状と同一と看做すことのできる形状に定められる。逆に、複数の磁気データの示す座標の分布形状から、楕円面の形状を特定することが困難な場合、複数の磁気データの示す座標と楕円面との誤差を最小化することにより生成される楕円面の形状は、誤差の表現形式に依存した形状となる。この発明によれば、複数の磁気データの示す座標と楕円面との誤差を第1軸成分の二乗値に基づいて表現したときに当該誤差を最小化する第1楕円面、複数の磁気データの示す座標と楕円面との誤差を第2軸成分の二乗値に基づいて表現したときに当該誤差を最小化する第2楕円面、及び、複数の磁気データの示す座標と楕円面との誤差を第3軸成分の二乗値に基づいて表現したときに当該誤差を最小化する第3楕円面のうち、少なくとも2つの楕円面を生成する。すなわち、初期楕円面生成部の生成する少なくとも2つの楕円面の各々は、異なる形式で表現された誤差を各々最小化することで定められる。従って、複数の磁気データの示す座標の分布形状から楕円面の形状を特定することが可能な場合、初期楕円面生成部の生成する少なくとも2つの楕円面は、複数の磁気データの示す座標の分布形状により特定される楕円面の形状と同一と看做すことのできる形状に定められる。一方、複数の磁気データの示す座標の分布形状から、楕円面の形状を特定することが困難な場合、初期楕円面生成部の生成する少なくとも2つ楕円面は、誤差の表現形式に依存した異なる形状を有する。この場合、少なくとも2つ楕円面の各々の中心点は、同一の座標と看做すことはできず、初期楕円面補正行列は生成されない。このように、本発明によれば、複数の磁気データの示す座標の分布の形状から、楕円面の形状を特定することが困難な場合に、不適切な初期楕円面補正行列が生成されることを防止することができる。また、上述した地磁気測定装置は、3次元の変数ベクトルの示す3軸の座標を起点として前記磁気データの示す3軸の座標を表す3次元のベクトルを第1変数ベクトルとし、前記第1変数ベクトルを、変数行列により変換して得られる3次元のベクトルを第2変数ベクトルとし、前記第2変数ベクトルの示す3軸の座標を表すデータを変換後データとし、複数の前記変換後データの各々が示す3軸の座標と、前記変数ベクトルの示す3軸の座標を中心とする球面との誤差を表し、前記変数行列の各成分と、前記変数ベクトルの各要素とを変数とする関数を、楕円面最適化関数としたとき、前記楕円面最適化関数の変数の初期値として、前記初期楕円面補正行列の各成分と、前記初期中心点の示す3軸の座標とを適用し、その後、前記楕円面最適化関数の変数の値を逐次更新することにより、前記楕円面最適化関数を最小化する解として、最適楕円面補正行列と、最適中心点の示す3軸の座標とを算出する、楕円面最適補正値生成部と、前記3次元磁気センサから出力される磁気データの示す3軸の座標を、前記最適中心点の示す3軸の座標から表した3次元のベクトルを、前記最適楕円面補正行列により変換することにより、地磁気の向きを算出する、地磁気算出部と、を更に備えることを特徴とすることが好ましい。この発明によれば、楕円面最適化関数の変数の初期値として、初期楕円面補正行列の各成分と、初期中心点の示す座標とを適用する。前述の通り、初期楕円面補正行列及び初期中心点は、初期楕円面生成部が生成する少なくとも2つ楕円面に基づいて生成される値であり、複数の磁気データの示す座標の分布形状を正確に表現する楕円面に基づいて定められる値である。一方、楕円面最適化関数を最小化する非線形最適化演算は、複数の磁気データの示す座標との誤差を最小化する楕円面を算出するための演算である。つまり、楕円面最適化関数の大局的最適解は、複数の磁気データの示す座標の分布形状を最も正確に表現する楕円面の形状を表す行列及び当該楕円面の中心点となる。従って、初期楕円面補正行列及び初期中心点は、楕円面最適化関数の大局的最適化に近い値を有する正確な値である。すなわち、本発明に係る地磁気測定装置は、大局的最適解に近い正確な値を、非線形最適化演算の初期値として適用することにより、非線形最適化演算が、局所最適解を算出することを防止し、正確な地磁気の向きを算出することができる。また、上述した地磁気測定装置において、前記3次元磁気センサから出力される磁気データの示す3軸の座標を、前記初期中心点の示す3軸の座標から表した3次元のベクトルを、前記初期楕円面補正行列により変換することにより、地磁気の向きを算出する、地磁気算出部を更に備えることを特徴とすることが好ましい。この発明によれば、簡易な計算により楕円補正を行うことができ、地磁気の向きの算出に係る計算負荷を軽減することが可能となる。<A.第1実施形態>以下、本発明の実施の形態を説明する。[1.3次元磁気センサが検出する磁界の概要]本実施形態では、3次元磁気センサが検出する磁界として、検出対象である地磁気の他に、3次元磁気センサが搭載される機器に備えられる部品が発する磁界である内部磁界、及び部品の備える軟磁性材料が機器外部からの磁界により磁化される結果生じる着磁磁界の存在を想定する。以下、図1乃至図5を用いて、本実施形態において想定するこれら3種類の磁界の概要、及び、これらの磁界を3次元磁気センサが検出した場合に3次元磁気センサから出力されるベクトルデータについて説明する。図1は、測定対象である地磁気Bg、3次元磁気センサの搭載された機器1が備える部品2が発する内部磁界Bi、及び、部品2のうち軟磁性材料21が発する着磁磁界Bmについて説明する図である。地磁気Bgは、磁極北に向かう一様な向き及び大きさを有する磁界である。地磁気Bgの向き及び大きさは、厳密には地域によって異なるが、例えば、異なる都市に移動する等の大きな移動をしない場合には、一様な向き及び大きさを有する。内部磁界Biは、機器1に備えられた部品2の発する磁界であり、機器1から見て一定の向き及び一定の大きさを有する磁界である。すなわち、機器1の姿勢がどのように変化しても、内部磁界Biは、3次元磁気センサ60によって、一定の向き及び大きさを有する磁界として検知される。着磁磁界Bmは、機器1の外部の物体から発せられる磁界すなわち、地磁気Bgの影響により軟磁性材料21が磁化される結果、軟磁性材料21が発する磁界である。従って、着磁磁界Bmは、地磁気Bgの向き及び大きさや、軟磁性材料21の材質、大きさ、形状等に依存して、向き及び大きさを変化させる磁界である。説明の便宜上、図1に示すような地上座標系ΣG及びセンサ座標系ΣSを導入する。図1に記載された各ベクトルの左上に付された添字Gは、当該ベクトルが地上座標系ΣGにおいて表現されたベクトルであることを意味する。地上座標系ΣGは、地上に固定された座標系であり、地上の任意の一点を原点とし、互いに直交する3つの方向、例えば、東、北、及び鉛直上向きを、それぞれx軸、y軸、及びz軸とする座標系である。センサ座標系ΣSは3次元磁気センサ60に固定された座標系であり、3次元磁気センサ60の有する3つの各々のセンサモジュールの出力する値を、それぞれ、x軸第1軸、y軸第2軸、z軸第3軸上にプロットするように設けられた座標系である。すなわち、3次元磁気センサ60が出力する磁気データは、センサ座標系ΣSのベクトルデータとして表現される。なお、図1に示す姿勢μは、地上座標系ΣGにおけるセンサ座標系ΣSの各軸の向きすなわち、地上座標系ΣGにおける、3次元磁気センサ60の向きを表す。以下において、姿勢μを変化させた場合、地上座標系ΣG及びセンサ座標系ΣSの各々において、内部磁界Bi及び着磁磁界Bmの向きがどのように変化するかについて説明する。まず、図2及び図3を参照しながら、地上座標系ΣG及びセンサ座標系ΣSにおいて、内部磁界Bi及び地磁気Bgがどのように表されるかについて説明する。なお、図2及び図3においては、簡単のために、機器1は軟磁性材料21を備えず、着磁磁界Bmが存在しない場合を想定する。図2は、内部磁界Bi及び地磁気Bgの向き及び大きさを、地上座標系ΣGにおいて表した図である。機器1の姿勢μが、姿勢μ1から姿勢μ2に変化した場合、内部磁界GBiは、大きさは一定であるが、向きは姿勢μの変化に伴い変化する。一方、地磁気GBgの向き及び大きさは一定である。図3は、内部磁界Bi及び地磁気Bgの向き及び大きさを、センサ座標系ΣSにおいて表現した図である。具体的には、図3は、機器1の姿勢μをμ1~μNと変化させつつ磁界を測定した場合に、3次元磁気センサ60が出力するN個の磁気データq1~qNの示す座標を、センサ座標系ΣSにおいてプロットした図であるNは精度よくオフセットを導出するために必要な磁気データの規定測定回数を表す9以上の自然数。ここで、図3に記載された各ベクトルの左上に付された添字Sは、当該ベクトルがセンサ座標系ΣSにおいて表されたベクトルであることを意味する。センサ座標系ΣSにおいて、内部磁界Biは、一定の向き及び一定の大きさを持つベクトルSBiセンサ座標系ΣSの原点から中心点cOGへと向かうベクトルとして表現される。一方、地磁気Bgは、大きさは一定であるが、向きは3次元磁気センサ60の姿勢μに伴い変化する。つまり、地磁気Bgは、機器1の姿勢μに依存した向きの、一定の大きさのベクトルSBgμとして表現される。従って、ベクトルSBgμの起点を中心点cOGに配置して姿勢μを変化させた場合、ベクトルSBgμの終点は、中心点cOGを中心点とし、地磁気Bgの大きさを半径とする、球面SG上の座標を示す。センサ座標系ΣSにおいて、複数の磁気データq1~qNの示す座標の各々は、内部磁界SBi及び地磁気SBgの和を表すため、複数の磁気データq1~qNの示す座標は、球面SG上に分布する。なお、3次元磁気センサ60の測定値は、測定誤差を有するため、複数の磁気データq1~qNの示す座標は、厳密には、球面SGの近傍に確率的に分布する。従って、磁気データqiの示す座標から内部磁界SBiを引き算することにより、センサ座標系ΣSにおける地磁気SBgの向き及び大きさを算出することが可能となる。このように、検出対象である地磁気Bgの正確な向きを得るために、磁気データqiの示す座標から、3次元磁気センサ60の出力する内部磁界Biを表す球面SGの中心点cOGの示す座標を減算する処理を、補正処理と呼ぶ。また、補正処理において磁気データqiから取り除かれるベクトルをオフセットcOFFと呼ぶ。すなわち、オフセットcOFFは、内部磁界を表すベクトルSBiであり、センサ座標系ΣSにおいて、原点から球面SGの中心点cOGを示すベクトルとして表される。機器1が軟磁性材料21を備える場合、地磁気Bgの影響を受けて軟磁性材料21が磁化される結果、軟磁性材料21は着磁磁界Bmを発する。機器1の姿勢μを変化させた場合、着磁磁界Bmの向き及び大きさが、地上座標系ΣG及びセンサ座標系ΣSの各々においてどのように変化するかについて、図4及び図5を参照しながら説明する。図4は、地上座標系ΣGにおける、着磁磁界Bmの向き及び大きさ例示した図である。図4は、機器1が、センサ座標系ΣSのx軸に平行な長辺211a及び211bと、センサ座標系ΣSのy軸に平行な短辺212a及び212bとを有する直方体の軟磁性材料21を備え、軟磁性材料21が、センサ座標系ΣSのx軸上に位置するように配置されている場合を例示している。着磁磁界Bmは、軟磁性材料21が地磁気Bgにより磁化された結果として生じる磁界であり、機器1の姿勢μ及び軟磁性材料21の材質、大きさ、形状等に依存して、向き及び大きさを変化させる磁界である。機器1の姿勢μが姿勢μ1から姿勢μ2に変化した場合、着磁磁界GBmの向き及び大きさは、GBmμ1からGBmμ2へと変化する。例えば、機器1の姿勢μが姿勢μ1の場合、軟磁性材料21は、軟磁性材料21の短辺212aから他方の短辺212bへと向かう着磁磁界GBmμ1を発生させ、機器1の姿勢μが姿勢μ2の場合、軟磁性材料21は、軟磁性材料21の長辺211aから他方の長辺211bへと向かう着磁磁界GBmμ2を発生させる。3次元磁気センサ60が検出する着磁磁界GBmμの向き及び大きさは、姿勢μと、センサ座標系ΣSにおける軟磁性材料21の位置SPmとに依存する。例えば、図4の場合、3次元磁気センサ60は、着磁磁界GBmμ1を、地磁気GBgμ1と同じ方向を向いた磁界として測定する。また、3次元磁気センサ60は、着磁磁界GBmμ2を、地磁気GBgμ2とは逆の方向を向いた磁界として測定する。なお、軟磁性材料21は、地磁気Bgの他に、部品2が発するセンサ座標系ΣSから見て一定の向き及び大きさの磁界により磁化される。センサ座標系ΣSから見て一定の向き及び大きさの磁界により磁化される結果として軟磁性材料21が発する磁界は、機器1の姿勢μが変化した場合であっても一定の方向及び大きさを有する。このような、軟磁性材料21が磁化された結果発する磁界のうち、姿勢μが変化しても一定の向き及び大きさを有する磁界は、前述した内部磁界Biに包含されるものとする。図5は、機器1が姿勢μ1となるときに測定される磁気データq1と、姿勢μ2となるときに測定される磁気データq2とを、センサ座標系ΣSにおいてプロットした図である。磁気データq1は、中心点cOGの示す座標を起点とし、地磁気SBgμ1と同じ向きを有する着磁磁界SBmμ1と、地磁気SBgμ1とを加算したベクトルSBEμ1により示される座標である。よって、磁気データq1は球面SGの外側に存在する。一方、磁気データq2は、中心点cOGの示す座標を起点とし、地磁気SBgμ2と逆向きの着磁磁界SBmμ2と、地磁気SBgμ2とを加算したベクトルSBEμ2により示される座標である。よって、磁気データq2は球面SGの内側に存在する。すなわち、磁気データq1及びq2は、球面SGを、ベクトルSBgμ1方向に引き延ばし、且つ、ベクトルSBgμ2方向に縮めることで得られる楕円面VE上に分布する。このように、3次元磁気センサが、軟磁性材料を備える機器に組み込まれている場合、軟磁性材料が地磁気等の機器の外部からの磁界により磁化される結果として生じる着磁磁界の影響により、3次元磁気センサの測定する複数の磁気データの示す座標が球面近傍に分布せず、楕円面近傍に分布する。このような、軟磁性材料が磁化された結果生じる磁界の影響により、複数の磁気データの示す座標が楕円面近傍に分布する現象は、ソフトアイアン効果と呼ばれている。ソフトアイアン効果が生じる場合に、3次元磁気センサ60が出力する複数の磁気データq1~qNの示す座標について、図6を参照しつつ説明する。図6は、3次元磁気センサ60の姿勢μをμ1~μNNは精度よくオフセットを導出するために必要な磁気データの規定測定回数を表す9以上の自然数と変化させつつ磁界を測定した場合に、3次元磁気センサ60が出力するN個の磁気データq1~qNの示す座標を、センサ座標系ΣSにおいてプロットした図である。図6は、ソフトアイアン効果により、複数の磁気データq1~qNの示す座標が、中心点cOGを中心とする楕円面VE上に分布する場合を想定している。なお、図6では、3次元磁気センサ60の測定誤差を考慮していないが、測定誤差を考慮する場合、複数の磁気データq1~qNの示す座標は、楕円面VE上には分布せず、楕円面VEの近傍に確率的に分布する。つまり、楕円面VEは、複数の磁気データq1~qNの示す座標との誤差を最小化するようにして定められる楕円面である。楕円面VEの主軸を、長い順番にLE1、LE2、及びLE3とし、これら3つの主軸の長さをrE1、rE2、及びrE3とする但し、rE1≧rE2≧rE3>0とする。また、球面SGの半径をrGとする。このとき、中心点cOGから磁気データqiの表す座標を示すベクトルSBEμiは、地磁気を表すベクトルSBgμiのうち、主軸LE1に平行な成分をrE1/rG倍したベクトル、主軸LE2に平行な成分をrE2/rG倍したベクトル、及び主軸LE3に平行な成分をrE3/rG倍したベクトルの和を表すベクトルとなる。従って、中心点cOGから磁気データqiの表す座標を示すベクトルSBEμiの向きと、地磁気を表すベクトルSBgμiの向きとは異なる。また、ベクトルSBEμi及びベクトルSBEμjとのなす角つまり、2つの磁気データqi及びqjの示す座標を、中心点cOGから見たときのなす角と、地磁気を表すベクトルSBgμi及びベクトルSBgμjのなす角とは、異なる。このように、ソフトアイアン効果が生じている場合、磁気データqiの座標から中心点cOGの座標を減算しても、地磁気SBgμiの向きを正確に求めることはできない。そこで、本実施形態では、図7に示すように、楕円面VE上の座標を、半径1の球面SE上の座標へと変換する楕円面補正行列TEを算出し、楕円面補正行列TEによって、磁気データqiの示す座標を、変換後磁気データsiにより表される球面SE上の座標へと変換する。中心点cOGから変換後磁気データsiの表す座標を示すベクトルSBsμiは、地磁気を表すベクトルSBgμiと同じ方向を向く。従って、変換後磁気データsiの示す座標から、中心点cOGの示す座標を減算することにより、地磁気を表すベクトルSBgμiの向きを求めることができる。このように、地磁気Bgの向きを算出するために、楕円面近傍に分布する複数の磁気データの示す座標を、当該楕円面と同一の中心点を有する半径1の球面近傍に分布する複数の座標へと変換する処理を、楕円補正と呼ぶ。楕円面補正行列TEが行う、楕円面VE上の磁気データqiの示す座標から、球面SE上の変換後磁気データsiの示す座標への座標変換は、以下の式1に表される。ここで、楕円面補正行列TEは、以下の式2に示される3行3列の対称行列である。また、式3に示される3次元の変数ベクトルqは、磁気データqiの座標を表すための変数ベクトルであり、式4に示される3次元の変数ベクトルsは、変換後磁気データsiの座標を表すための変数ベクトルであり、式5に示される3次元の変数ベクトルcは、中心点cOGの座標すなわち、オフセットcOFFを表すための変数ベクトルである。なお、式1において、ベクトルq-cは、楕円面VEの中心点cOGをセンサ座標系ΣSの原点に平行移動させて得られる楕円面上の座標を示し、ベクトルs-cは、座標系ΣSの原点を中心とする半径1の球面上の座標を示す。前述の通り、楕円面補正行列TEは、楕円面VEの中心点cOGを原点とする座標系において、楕円面VE上の座標を、半径1の球面SE上の座標へと変換するための行列である。つまり、楕円面補正行列TEは、互いに直交する3つの固有ベクトルの各々と、楕円面VEの3つの主軸の各々とが平行になり、3つの固有ベクトルに対応する3つの固有値の各々と、楕円面VEの3つの主軸の各々の長さの逆数とが等しくなるように定められる。ここで、楕円面補正行列TEの有する3つの固有ベクトルをuT1、uT2、及びuT3とし、対応する固有値をλT1、λT2、及びλT3但し、λT1≧λT2≧λT3>0とする。このとき、固有ベクトルuT1は、主軸LE1と平行となり、固有ベクトルuT2は、主軸LE2と平行となり、固有ベクトルuT3は、主軸LE3と平行となるように定められる。また、固有値λT1は、主軸LE1の長さrE1の逆数に等しくなり、固有値λT2は、主軸LE2の長さrE2の逆数に等しくなり、固有値λT3は、主軸LE3の長さrE3の逆数に等しくなるように定められる。すなわち、楕円面補正行列TEは、任意のベクトルに対して、固有ベクトルuT1方向の成分を固有値λT1倍に伸縮し、固有ベクトルuT2方向の成分を固有値λT2倍に伸縮し、固有ベクトルuT3方向の成分を固有値λT3倍に伸縮する行列である。なお、楕円面補正行列TEの有する3つの固有値をλT1、λT2、及びλT3は、いずれも正の値であり、楕円面補正行列TEは、正定値行列である。ところで、図7に示すように、中心点cOGから変換後磁気データsiにより表される座標を示すベクトルSBsμiと、地磁気を表すベクトルSBgμiとの間には、ずれ角ψが生じる場合がある。この場合、変換後磁気データsiから、ベクトルSBgμiの正確な向きを算出することはできない。しかし、ずれ角ψは、軟磁性材料21と3次元磁気センサ60との相互位置関係すなわち、ベクトルSPmの向き及び大きさに依存する値である。従って、ベクトルSPmよりずれ角ψを特定することが可能であり、特定されたずれ角ψと、複数の変換後磁気データsiとから、正確な地磁気Bgの向きを算出することが可能である。また、軟磁性材料21の配置を工夫することにより、ずれ角ψを小さく抑えることも可能である。以下において、楕円面VEの形状を求め、楕円面補正行列TEを算出することにより楕円補正を行い、正確な地磁気Bgの方向を算出する方法について述べる。[2.機器構成及びソフトウェア構成]図8は、本発明の第1実施形態に係る機器1の構成を示すブロック図である。機器1は、各種の構成要素とバスを介して接続され装置全体を制御するCPU10、CPU10の作業領域として機能するRAM20、各種プログラムやデータを記憶したROM30、通信を実行する通信部40、画像を表示する表示部50、及び磁気を検出して磁気データを出力する3次元磁気センサ60を備える。3次元磁気センサ60は、X軸地磁気センサ61、Y軸地磁気センサ62、及びZ軸地磁気センサ63を備える。各センサは、MI素子磁気インピーダンス素子、MR素子磁気抵抗効果素子などを用いて構成することができる。地磁気センサI/F64は、各センサから出力信号をAD変換して磁気データqを出力する。この磁気データqは、X軸地磁気センサ61、Y軸地磁気センサ62、及びZ軸地磁気センサ63からの出力の各々を、センサ座標系ΣSの、x軸、y軸及びz軸の3成分によって示す、センサ座標系ΣS上のベクトルデータである。CPU10、RAM20、3次元磁気センサ60、及び磁気データ処理プログラム70は、3次元磁気センサ60が検出し出力する磁気データqに基づき、正確な地磁気の向きを示す地磁気データを算出する地磁気測定装置として機能する。表示部50は、CPU10が磁気データ処理プログラム70を実行することにより算出した地磁気の向きを、方位情報として矢印等により表示する。なお、磁気データ処理プログラム70は、地図アプリケーション等との連携を想定するものであってもよく、表示部50は、地磁気の向きを示す方位情報である矢印等を地図上に表示してもよい。図9は、地磁気測定装置のうち、CPU10が磁気データ処理プログラム70を実行することにより実現される機能を表す、機能ブロック図である。地磁気測定装置は、複数の磁気データq1~qNを蓄積する蓄積部100と、最適中心点cEOPの座標及び最適楕円面補正行列TOPを算出する楕円面補正部200と、磁気データqi、最適中心点cEOP、及び最適楕円面補正行列TOPに基づいて地磁気Bgの方向を算出する地磁気算出部600と、を備える。ここで、最適中心点cEOPとは、複数の磁気データq1~qNの示す座標との誤差を最小化するように定められる楕円面である最適楕円面VEOPの中心点である。また、最適楕円面補正行列TOPとは、最適楕円面VEOP上の座標を、最適中心点cEOPを中心とする球面SEOP上の座標へと変換するための、3行3列の対称行列である。蓄積部100は、3次元磁気センサ60から順次出力される磁気データq1~qNをバッファBU1に蓄積するNは、精度よくオフセットを導出するために必要な磁気データの規定測定回数を表す9以上の自然数。バッファBU1はRAM20によって形成される。楕円面補正部200は、楕円面初期補正値生成部300と、楕円面最適補正値生成部400とを備える。楕円面初期補正値生成部300は、蓄積部100に蓄積された複数の磁気データq1~qNに基づいて初期楕円面補正行列T0及び初期中心点cE0の座標を算出する。ここで、初期中心点cE0とは、蓄積部100に蓄積された複数の磁気データq1~qNの示す座標を近傍に有する初期楕円面VE0の中心点である。また、初期楕円面補正行列T0とは、初期楕円面VE0上の座標を、初期中心点cE0を中心とする球面SE0上の座標へと変換するための3行3列の対称行列である。楕円面最適補正値生成部400は、楕円面初期補正値生成部300が出力する初期楕円面補正行列T0及び初期中心点cE0の座標に基づいて、複数の磁気データq1~qNの各々が示す座標と楕円面との誤差を最小化する最適楕円面VEOPの中心点である最適中心点cEOPの座標と、最適楕円面VEOP上の座標から最適中心点cEOPを中心点に有する球面SEOP上の座標への座標変換を表す最適楕円面補正行列TOPとを算出する。複数の磁気データq1~qNの示す座標と最適楕円面VEOPとの誤差が最小化され0となる場合、楕円面VEと最適楕円面VEOPとは一致し、最適中心点cEOPと中心点cOGすなわち、内部磁界Biの示す座標とは一致する。地磁気算出部600は、3次元磁気センサ60から出力される磁気データqiの示す座標に対して、楕円面補正行列TEとオフセットcOFFとを用いた楕円補正を施し、センサ座標系ΣSにおける地磁気SBgの向き厳密には、ベクトルSBsμiの向きを算出する。具体的には、地磁気算出部600は、オフセット採用部610と、地磁気ベクトル計算部620とを備える。オフセット採用部610は、楕円面補正行列TEとして最適楕円面補正行列TOPを採用するとともに、オフセットcOFFとして最適中心点cEOPの座標を示すベクトルを採用する。また、地磁気ベクトル計算部620は、3次元磁気センサ60から出力される磁気データqiに対して、楕円面補正行列TE及びオフセットcOFFを用いた楕円補正を行い、地磁気SBgの向きを算出する。以下において、楕円面初期補正値生成部300、楕円面最適補正値生成部400、及び、地磁気算出部600の詳細について、述べる。[3.初期楕円面の生成]図10は、楕円面初期補正値生成部300の機能構成を示す機能ブロック図である。本実施形態では、複数の磁気データq1~qNに基づいて初期楕円面VE0を算出する際に、複数の磁気データq1~qNの示す座標を近傍に有する第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzを生成し、これら3つの楕円面に基づいて初期楕円面VE0を生成する。以下、本実施形態における初期楕円面VE0の生成方法について具体的に説明する。楕円面初期補正値生成部300は、第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzの各々の係数行列Dxx、Dyy、Dzz及び中心点cxx、cyy、czzの座標を算出する初期楕円面生成部310、第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzの各々の係数行列及び中心点に基づいて初期楕円面補正行列T0及び初期中心点cE0の座標を算出することが適切であるか否かを判断する初期楕円面判定部320、及び、第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzの各々の係数行列及び中心点の座標に基づいて、初期楕円面補正行列T0及び初期中心点cE0の座標を算出する初期補正値生成部330を備える。初期楕円面生成部310は、蓄積部100に蓄積された複数の磁気データq1~qNに基づいて、第1楕円面Vxxの形状を表す第1楕円面係数行列Dxxと、第1楕円面Vxxの中心点cxxの座標とを算出する第1楕円面生成部311、複数の磁気データq1~qNに基づいて、第2楕円面Vyyの形状を表す第2楕円面係数行列Dyyと、第2楕円面Vyyの中心点cyyの座標とを算出する第2楕円面生成部312、及び、複数の磁気データq1~qNに基づいて、第3楕円面Vzzの形状を表す第3楕円面係数行列Dzzと、第3楕円面Vzzの中心点czzの座標とを算出する第3楕円面生成部313、を備える。以下において、第1楕円面係数行列Dxx、第2楕円面係数行列Dyy、及び第3楕円面係数行列Dzzと、中心点cxxの座標、中心点cyyの座標、及び中心点czzの座標の算出方法について述べる。3次元磁気センサ60より出力される磁気データqの示す座標を表す変数を、式3で表す場合、磁気データqを表面上に有する楕円面の方程式楕円面方程式は、以下の式6で表される。なお、式6は楕円面を表すので、式6に現れる係数θxx、θyy、及びθzzはいずれも正の値となる。式6に示される楕円面方程式は、以下の式7に変形される。複数の磁気データq1~qNの示す座標の全てが、式6に示された楕円面上に位置する場合、以下の式8が成立する。但し、ベクトルθXXは、式9に示すように、式7の9つの係数を並べた9次元のベクトルである。また、行列Rxxは、式10に示すように、式13に示す9次元のベクトルQxxに対して、式11に示す複数の磁気データq1~qNの示す座標の各々を代入して得られるN個のベクトルを転置したうえで各行に並べたN行9列の行列である。ベクトルWxxは、式12に示すように、複数の磁気データq1~qNの示す座標のうち、x成分の2乗値にマイナスを掛けた値を各成分に有する、N次元のベクトルである。式8は、ベクトルθXXの各要素を変数とする連立一次方程式である。従って、式8をベクトルθXXについて解くことで、式7の係数が決定され、複数の磁気データq1~qNの示す座標を表面上に有する楕円面方程式を特定できる。しかし、3次元磁気センサ60の測定誤差等を考慮すると、複数の磁気データq1~qNの示す各々の座標の全てが、式7に示す楕円面上に正確に一致する位置に存在することはない。従って、式8は解を有さず、式8の解としてベクトルθXXを算出することはできない。そこで、本実施形態では、式8の解として尤もらしいベクトルθXXを、統計的手法を用いることにより算出する。例えば、式7の右辺に現れる8つの項xy、xz、y2、yz、z2、x、y、及びzをそれぞれ独立変数と看做し、式7の左辺に現れるx2を従属変数と看做したうえで、最小二乗法を用いて式14に示す正規方程式を導出し、その解として、ベクトルθXXを求める。この正規方程式の解として表されるベクトルθXXは、行列RxxTRxxが正則なときに、式15で表すことができる。式15により算出されたベクトルθXXを係数として式7に適用することにより表される楕円面を、第1楕円面Vxxと呼ぶ。ここで、図11Aに示すように、xy、xz、y2、yz、z2、x、y、及びzの各々を変数として表す8本の軸を有する8次元の空間πxxと、x2の値を変数として表す第1評価軸ξ1とからなる、9次元の空間Ω上に、複数の磁気データq1~qNの各々をプロットする。このとき、楕円面Vxxは、楕円面Vxxと、複数の磁気データq1~qNが空間Ω上で示す座標との、第1評価軸ξ1方向の誤差を最小化する立体空間Ω上の8次元平面として求められる。すなわち、楕円面Vxxは、複数の磁気データq1~qNを、空間πxx上にプロットして得られる複数の8次元のベクトルqπxx1~qπxxNの各要素を式7の右辺に代入して得られる値qxx1~qxxNと、複数の磁気データq1~qNの示す座標のうちx軸成分の2乗値x12~xN2との誤差を最小化する立体として定められる。式7に示される第1楕円面Vxxの方程式は、式17に示す第1楕円面係数行列Dxxを用いて、式16に変形される。このとき、第1楕円面Vxxの中心点cxxの座標は、式18で表される。第1楕円面生成部311は、以上に示したように、第1楕円面係数行列Dxx及び第1楕円面Vxxの中心点cxxを算出し、これらを出力する。なお、式16が楕円面を表すための必要条件は、第1楕円面係数行列Dxxが正定値となることである。次に、式6に示される楕円面方程式を、以下の式19に変形する。式19に示される方程式は、ベクトルθYYの各要素を変数とする連立一次方程式である式20に変形される。式8と同様に、式20の解としてベクトルθYYを算出することはできないため、式20の解として尤もらしい値として、ベクトルθYYを算出する。具体的には、式24に示す正規方程式において、行列RyyTRyyが正則なときには、ベクトルθYYは式25により算出される。式25により特定されたベクトルθYYを係数として式19に適用することにより表される楕円面を、第2楕円面Vyyと呼ぶ。なお、ベクトルθYYは、式21に示す9次元のベクトルであり、行列Ryyは、式22に示すN行9列の行列であり、ベクトルWyyは、式23に示すN次元のベクトルである。ここで、図11Bに示すように、x2、xy、xz、yz、z2、x、y、及びzの各々を変数として表す8本の軸を有する8次元の空間πyyと、y2の値を変数として表す第2評価軸ξ2とからなる、9次元の空間Ω上に、複数の磁気データq1~qNの各々をプロットする。このとき、楕円面Vyyは、楕円面Vyyと、複数の磁気データq1~qNが空間Ω上で示す座標との、第2評価軸ξ2方向の誤差を最小化する立体空間Ω上の8次元平面として求められる。すなわち、楕円面Vyyは、複数の磁気データq1~qNを、空間πyy上にプロットして得られる複数の8次元のベクトルqπyy1~qπyyNの各要素を式19の右辺に代入して得られる値qyy1~qyyNと、複数の磁気データq1~qNの示す座標のうちy軸成分の2乗値y12~yN2との誤差を最小化する立体として定められる。式19に示される第2楕円面Vyyの方程式は、式27に示す第2楕円面係数行列Dyyを用いて、式26に変形される。このとき、第2楕円面Vyyの中心点cyyの座標は、式28で表される。第2楕円面生成部312は、以上に示したように、第2楕円面係数行列Dyy及び第2楕円面Vyyの中心点cyyを算出し、これらを出力する。なお、式26が楕円面を表すための必要条件は、第2楕円面係数行列Dyyが正定値となることである。次に、式6に示される楕円面方程式を、以下の式29に変形する。式29に示される方程式は、ベクトルθZZの各要素を変数とする連立一次方程式である式30に変形される。式8と同様に、式29の解としてベクトルθZZを算出することはできないため、式29の解として尤もらしい値として、ベクトルθZZを算出する。具体的には、式34に示す正規方程式において、行列RzzTRzzが正則なときには、ベクトルθZZは式35により算出される。式35により特定されたベクトルθZZを係数として式29に適用することにより表される楕円面を、第3楕円面Vzzと呼ぶ。なお、ベクトルθZZは、式31に示す9次元のベクトルであり、行列Rzzは、式32に示すN行9列の行列であり、ベクトルWzzは、式33に示すN次元のベクトルである。ここで、図11Cに示すように、x2、xy、xz、y2、yz、x、y、及びzの各々を変数として表す8本の軸を有する8次元の空間πzzと、z2の値を変数として表す第3評価軸ξ3とからなる、9次元の空間Ω上に、複数の磁気データq1~qNの各々をプロットする。このとき、楕円面Vzzは、楕円面Vzzと、複数の磁気データq1~qNが空間Ω上で示す座標との、第3評価軸ξ3方向の誤差を最小化する立体空間Ω上の8次元平面として求められる。すなわち、楕円面Vzzは、複数の磁気データq1~qNを、空間πzz上にプロットして得られる複数の8次元のベクトルqπzz1~qπzzNの各要素を式29の右辺に代入して得られる値qzz1~qzzNと、複数の磁気データq1~qNの示す座標のうちz軸成分の2乗値z12~zN2との誤差を最小化する立体として定められる。式29に示される第3楕円面Vzzの方程式は、式37に示す第3楕円面係数行列Dzzを用いて、式36に変形される。このとき、第3楕円面Vzzの中心点czzの座標は、式38で表される。第3楕円面生成部313は、以上に示したように、第3楕円面係数行列Dzz及び第3楕円面Vzzの中心点czzを算出し、これらを出力する。なお、式36が楕円面を表すための必要条件は、第3楕円面係数行列Dzzが正定値となることである。このように、初期楕円面生成部310は、第1楕円面係数行列Dxx、第2楕円面係数行列Dyy、第3楕円面係数行列Dzz、中心点cxxの座標、中心点cyyの座標、及び中心点czzの座標を算出し、これらを出力する。初期楕円面判定部320は、図10に示すように、初期楕円面係数行列判定部321と、初期楕円面中心点判定部322とを備える。初期楕円面係数行列判定部321及び初期楕円面中心点判定部322には、第1楕円面係数行列Dxx、第2楕円面係数行列Dyy、第3楕円面係数行列Dzz、中心点cxxの座標、中心点cyyの座標、及び中心点czzの座標が供給される。初期楕円面係数行列判定部321は、第1楕円面係数行列Dxx、第2楕円面係数行列Dyy、及び第3楕円面係数行列Dzzの全てが正定値であるという条件第1の条件が満たされるか否かを判定する。また、初期楕円面中心点判定部322は、式39に示すように中心点cxxと中心点cyyとの距離が第1閾値Δc以下であり、且つ、式40に示すように中心点cyyと中心点czzとの距離が第1閾値Δc以下であり、且つ、式41に示すように中心点czzと中心点cxxとの距離が第1閾値Δc以下であるという条件第2の条件が満たされるか否かを判定する。第1の条件に係る判定結果が肯定であり、且つ、第2の条件に係る判定結果が肯定である場合、初期楕円面判定部320は、第1楕円面係数行列Dxx、第2楕円面係数行列Dyy、第3楕円面係数行列Dzz、中心点cxxの座標、中心点cyyの座標、及び中心点czzの座標を出力する。一方、第1の条件または第2の条件の判定結果が否定である場合、地磁気測定装置は、処理を中断する。なお、本実施形態では、初期楕円面判定部320は、初期楕円面係数行列判定部321及び初期楕円面中心点判定部322を備え、第1の条件及び第2の条件の双方を満たすか否かを判定するが、本発明はこのような判定方法に限定されるものではない。例えば、初期楕円面判定部320は、初期楕円面係数行列判定部321を備えずに構成されてもよい。この場合、初期楕円面判定部320は、第1の条件についての判定を行わず、第2の条件についての判定のみを行い、判定結果が肯定である場合に、第1楕円面係数行列Dxx、第2楕円面係数行列Dyy、第3楕円面係数行列Dzz、中心点cxxの座標、中心点cyyの座標、及び中心点czzの座標を出力してもよい。初期補正値生成部330は、初期楕円面判定部320が行う判定の判定結果が肯定である場合、初期楕円面判定部320からの出力に基づいて、初期楕円面補正行列T0及び初期中心点cE0を算出する。ここで、初期楕円面補正行列T0とは、図12に示すように、初期中心点cE0を中心点とする初期楕円面VE0上の座標を、初期中心点cE0を中心点とする球面SE0上の座標へと変換する行列である。初期楕円面VE0は、第1楕円面Vxx、第2楕円面Vyy、及び、第3楕円面Vzzのうち、少なくとも1つの楕円面に基づいて定められる楕円面であり、複数の磁気データq1~qNの示す座標を表面近傍に有する楕円面である。初期楕円面補正行列T0の具体的な算出方法を説明するにあたり、まず、楕円面VEの形状が既知であると仮定したときの、楕円面補正行列TEの算出方法について説明する段落0043及び図7参照。楕円面補正行列TEは、楕円面VE上の座標を、楕円面VEの中心点cOGを中心とする半径1の球面SE上の座標に変換する行列であり、楕円面VEの形状を表す楕円面係数行列Dと、楕円面VEの中心点cOGとに基づき算出される。式6に示した楕円面VEを表す楕円面方程式は、以下の式43に示す楕円面係数行列Dを用いて、以下の式42に変形できる。また、楕円面VEの中心点cOGの示す座標は、以下の式44により表される。ここで、M行M列の正定値対称行列Gと、M行M列の正定値対称行列Hとの間に、式45の関係が成立する場合、行列Gを行列Hの平方根行列と呼ぶ。以下では、行列Hの平方根行列Gを、式46のように、行列の2分の1乗として表現する。このとき、行列Hの平方根行列Gは、式47により求められる。但し、行列U及び行列Λは、式48に示すように、行列Hを対角化することで算出される行列である。具体的には、行列Λは、式49に示すように、行列HのM個の正の固有値λH1~λHMを対角成分とするM行M列の行列であり、行列Uは、行列Hの固有値λH1~λHMに対応する固有ベクトルを正規化したうえで各列に並べたM行M列の回転行列である。楕円面係数行列Dと、楕円面補正行列TEとの関係は、以上で定義した平方根行列を用いて、以下の式50で表される。なお、値rDは、以下の式51に示すように、複数の磁気データq1~qNの示す座標を楕円面係数行列Dの平方根行列により変換した複数の座標と、中心点cOGとの距離の平均値を表す値である。楕円面係数行列Dの固有値λD1、λD2、及びλD3の各々は、楕円面VEの主軸の長さrE1、rE2、及びrE3の各々の2乗値の逆数に等しい。よって、楕円面補正行列TEの固有値λT1、λT2、及びλT3の各々は、楕円面VEの主軸の長さrE1、rE2、及びrE3の各々の逆数に等しくなる。従って、楕円面補正行列TEにより、楕円面VE上の座標を、半径1の球面SE上の座標へと変換することができる。次に、初期楕円面補正行列T0、及び、初期中心点cE0の座標の算出方法について説明する。初期楕円面補正行列T0は、第1楕円面係数行列Dxx、第2楕円面係数行列Dyy、及び第3楕円面係数行列Dzzのうち、少なくとも1つに基づいて算出される。また、初期中心点cE0の座標は、中心点cxx、中心点cyy、及び中心点czzのうち、少なくとも1つに基づいて算出される。前述した第2の条件を満たす場合、中心点cxx、中心点cyy、及び中心点czzのうち、任意の2点間の距離は、第1閾値Δcよりも近い。従って、第1閾値Δcの値を十分に小さくした場合、第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzは、いずれも複数の磁気データq1~qNの示す座標を近傍に有し、いずれの中心点も同一の座標と看做すことができるため、これら3つの楕円面は厳密には異なる形状を有するものの実質的に同一の楕円面と看做すことができる。この場合、初期楕円面VE0として、第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzのいずれを採用してもよい。本実施形態では、初期楕円面VE0として、第1楕円面Vxxを採用する。このとき、初期楕円面補正行列T0及び、初期中心点cE0の座標は、以下に示す式52及び式53により表される。初期補正値生成部330は、以上のように、初期楕円面補正行列T0と初期中心点cE0の座標とを生成し、これらを出力する。なお、本実施形態では、楕円面初期補正値生成部300は、蓄積部100から取得した複数の磁気データq1~qNを、楕円面最適補正値生成部400に対して出力するが、楕円面最適補正値生成部400は、蓄積部100より直接、複数の磁気データq1~qNを取得してもよい。また、本実施形態では、初期楕円面VE0として、第1楕円面Vxxを採用するが、本発明はこのような形態に限定するものではない。例えば、初期楕円面VE0として、第2楕円面Vyyを採用してもよい。このとき、初期楕円面補正行列T0は式54により表され、初期中心点cE0には中心点cyyが採用される。また、初期楕円面VE0として第3楕円面Vzzを採用してもよい。この場合、初期楕円面補正行列T0は式55により表され、初期中心点cE0には中心点czzが採用される。また、初期楕円面補正行列T0は、以下に示す式56により算出してもよい。この場合、初期中心点cE0の座標を、以下に示す式57または式58により算出してもよい。ところで、初期楕円面補正行列T0を算出する方法には、以下に示す方法以下、対比例と称するも存在する非特許文献2参照。具体的には、まず、式6に示す楕円面方程式を、x2の項、y2の項、またはz2の項の何れか1つの項で割り算をすることにより、楕円面方程式と等価な、以下の式59に示す連立一次方程式に変形する。次に、最小二乗法を用いて、式59から、以下の式60に示す正規方程式を算出する。行列RTRが正則なときには、楕円面の形状を示すベクトルθは、式60に示す正規方程式の解として、以下の式61により算出される。式61により算出されるベクトルθと、式43、式44、及び50とを用いて、初期楕円面補正行列T0及び初期中心点cE0を算出する。なお、例えば、式6に示す楕円面方程式を、z2の項により割り算することで、式59に示す連立一次方程式を算出する場合、ベクトルθは、以下の式62に示す9次元のベクトルであり、行列Rは、以下の式63に示す9次元のベクトルQに対して式11に示す複数の磁気データq1~qNの示す座標を代入して得たベクトルを転置して各行に並べて生成される以下の式64に示すN行9列の行列であり、ベクトルWは、以下の式65に示す9次元のベクトルである。対比例は、楕円面方程式を、x2の項、y2の項、またはz2の項の何れか1つの項で割り算をすることで生成される連立一次方程式により、楕円面の形状を決定するものである。つまり、対比例は、第1評価軸ξ1、第2評価軸ξ2、及び第3評価軸ξ3のうち、1つの評価軸のみを用いて、1つの楕円面を算出することと等しい。この場合、複数の磁気データq1~qNの示す座標の分布形状が同一であっても、選択される評価軸が相違すれば、算出される楕円面の形状は異なる。しかし、対比例は、1つの楕円面のみを算出するため、仮に2つ以上の評価軸を用いて2つ以上の楕円面を算出したならば把握できるはずの複数の楕円面の形状の相違について、把握することができない。例えば、2つ以上の楕円面が算出されることを前提とする、第2の条件を用いた判定を行うことができない。従って、図13に示す例のように、本実施形態に係る方法により、第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzのうち、少なくとも2つの楕円面を算出していたならば、第2の条件についての判定結果が否定となるような場合であっても、対比例によれば初期楕円面補正行列T0が算出されることになる。このように、対比例は、複数の磁気データq1~qNの示す座標の分布の形状から楕円面の形状を特定することが困難な場合であっても、複数の磁気データq1~qNの示す座標の分布形状を正確に表さない不適切な初期楕円面VE0に基づいた初期楕円面補正行列T0を生成する。これに対して、楕円面初期補正値生成部300は、初期楕円面補正行列T0及び初期中心点cE0を生成するにあたり、第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzの3つの楕円面を生成する。第1楕円面Vxxは、空間Ω上の第1評価軸ξ1を基準として複数の磁気データq1~qNの示す座標との誤差を最小化する楕円面であり、第2楕円面Vyyは、空間Ω上の第2評価軸ξ2を基準として複数の磁気データq1~qNの示す座標との誤差を最小化する楕円面であり、第3楕円面Vzzは、空間Ω上の第3評価軸ξ3を基準として複数の磁気データq1~qNの示す座標との誤差を最小化する楕円面である。すなわち、楕円面初期補正値生成部300は、空間Ω上の3つの異なる評価軸を用いて、3つの異なる楕円面を生成する。さらに、楕円面初期補正値生成部300は、第1の条件及び第2の条件を用いて、算出された3つの異なる楕円面が、互いに近い形状を有することを判定する。つまり、複数の磁気データq1~qNの示す座標の分布の形状が楕円面の形状とは大きく異なる形状であり、第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzのうち1つ以上が楕円面とは異なる形状となる場合、第1の条件は満たされない。また、図13に示す例のように、算出された3つの楕円面各々の中心点の相互間距離が離れている場合には、第2の条件は満たされない。このように、楕円面初期補正値生成部300は、初期楕円面補正行列T0を生成する際に、3つの異なる楕円面を生成したうえで、第1の条件及び第2の条件が満たされているか否かを判定するため、複数の磁気データq1~qNの示す座標の分布の形状から、楕円面の形状を特定することが困難な場合に、複数の磁気データq1~qNの示す座標の分布の形状とは異なる不正確な初期楕円面VE0に基づいた、不適切な初期楕円面補正行列T0が生成されることを防止する。[4.最適楕円面の生成]初期楕円面補正行列T0は、式66に示すように、初期中心点cE0を起点として磁気データqiの示す座標を終点とするベクトルqi-cEOを、初期楕円面VE0の3つの主軸方向に伸縮させることで、初期中心点cE0を起点として変換後磁気データs0iの示す座標を終点とするベクトルs0i-cEOに変換する行列である。磁気データqiの示す座標が、初期楕円面VE0上に存在する場合には、変換後磁気データs0iの示す座標は、初期中心点cE0を中心とする球面SE0上に位置する。図12に示すように、初期楕円面VE0は、複数の磁気データq1~qNの示す座標を近傍に有するように定められる楕円面であるが、複数の磁気データq1~qNの示す座標との誤差を最小化するように定められた楕円面ではない。従って、複数の磁気データq1~qNの示す座標と、初期楕円面VE0との間の誤差が大きい場合、初期楕円面VE0の中心点である初期中心点cE0の示す座標と、地磁気Bgを表す球面SGの中心点cOGの示す座標との間の誤差も大きくなる。この場合、磁気データqiの示す座標を初期楕円面補正行列T0及び初期中心点cE0により楕円補正した変換後磁気データs0iと、初期中心点cE0とに基づいて、地磁気Bgの正確な向きを算出することはできない。そこで、本実施形態では、複数の磁気データq1~qNの示す座標との誤差を最小化する最適楕円面VEOPに基づいて定められる最適楕円面補正行列TOPと、最適楕円面VEOPの中心点である最適中心点cEOPとに基づいて、磁気データqiの示す座標を楕円補正する。具体的には、本実施形態に係る地磁気測定装置は、楕円面補正行列TEとして最適楕円面補正行列TOPを採用するとともに、オフセットcOFFとして最適中心点cEOPの示す座標を採用することで、式1を、以下の式67に変形し、地磁気SBgの向きを表すベクトルsi-cEOPを算出する。以下、式67の右辺に現れるベクトルqi-cEOPを第1磁気ベクトルと称し、式67の左辺に現れるベクトルsi-cEOPを第2磁気ベクトルと称する。最適楕円面VEOPは、複数の磁気データq1~qNの示す座標との誤差を最小化する楕円面であるため、楕円面VEと、最適楕円面VEOPとは、同一の図形と看做すことができる。よって、最適中心点cEOPの示す座標と中心点cOGの示す座標オフセットcOFFとの誤差は、初期中心点cE0の示す座標と中心点cOGの示す座標との誤差に比べて小さく、最適中心点cEOPと中心点cOGとは、同一の座標を示すと看做すことができる。このような、最適楕円面補正行列TOPと、最適中心点cEOPとに基づいて、磁気データqiの示す座標を楕円補正することにより、正確な地磁気Bgの向きを得ることができる。楕円面最適補正値生成部400は、以下の式68に示す楕円面最適化関数fELT,cの値が最小化されるように、変数行列Tの各成分及び変数ベクトルcの各要素を逐次更新する非線形最適化演算を行い、楕円面最適化関数fELT,cの値が最小化されたときの変数行列T及び変数ベクトルcを、最適楕円面補正行列TOP及び最適中心点cEOPとして算出する。ここで、楕円面最適化関数fELT,cは、以下の式69に示す3行3列の対称行列である変数行列Tの各成分と、式5に示した変数ベクトルcの各要素とを変数とする関数であり、以下の式70として表すことができる。変数行列T及び変数ベクトルcの初期値には、初期楕円面補正行列T0及び初期中心点cE0が適用される。楕円面最適化関数fELT,cは、式68に示すように、変数ベクトルcの示す座標を起点として複数の磁気データq1~qNの示す座標を表す複数の第1変数ベクトルqi-cの各々を、変数行列Tにより変換することで得られる複数の第2変数ベクトルTqi-cの長さの平均値が、1からどの程度の差異を有するかを表す関数である。つまり、楕円面最適化関数fELT,cは、複数の第2変数ベクトルTqi-cを、変数ベクトルcの示す座標が起点となるように配置した場合、複数の第2変数ベクトルの各々が示す座標と、変数ベクトルcの示す座標を中心とする半径1の球面との誤差を表す。このとき、複数の第2変数ベクトルTqi-cの示す複数の座標を表すデータを、複数の変換後データsX1~sXNと呼ぶ。楕円面最適化関数fELT,cの値を最小化することにより、複数の変換後データsX1~sXNの示す座標と、変数ベクトルcの示す座標を中心とする半径1の球面との誤差を最小化することができ、このときの複数の変換後データsX1~sXNは、複数の変換後磁気データs1~sNを表す。なお、本実施形態では、説明の便宜上、第2変数ベクトルTqi-cを、変数ベクトルcの示す座標を起点として配置しているが、第2変数ベクトルTqi-cを、センサ座標系ΣSの原点が起点となるように配置しても構わない。すなわち、式68は、第2変数ベクトルTqi-cをセンサ座標系ΣSの原点が起点となるように配置した場合、第2変数ベクトルTqi-cの示す座標と、センサ座標系ΣSの原点を中心とする半径1の球面との誤差を表す。また、この場合、複数の変換後磁気データs1~sNは、原点を中心とする半径1の球面近傍に分布する。楕円面最適化関数fELT,cの示す値を最小化することにより最適楕円面補正行列TOP及び最適中心点cEOPを算出する非線形最適化演算は、公知の方法を適宜適用すればよい。例えば、非線形最適化演算として、ニュートン法を適用してもよい。ニュートン法、最急降下法等の、非線形最適化演算は、非線形関数の示す値を最適化最小化または最大化するように、非線形関数の変数の値を逐次更新する。そして、非線形最適化演算は、非線形関数の値及び変数の値が予め規定された停止規準を満たしたときに、変数の値の更新を停止し、このときの変数の値を最適解として採用する。なお、非線形最適化演算の停止規準としては、公知の基準を適宜適用すればよい。例えば、アルミホの規準Armijo’sruleを適用してもよい。非線形最適化演算は、非線形関数の最適解、すなわち非線形関数を最小化または最大化する大局的最適解を算出するための演算である。しかし、非線形最適化演算に適用される初期値が大局的最適解とは大きく異なる値である場合、非線形最適化演算の算出する最適解が、大局的最適解とは異なる局所的最適解となることがある。大局的最適解と大きく異なる初期値を適用した場合、大局的最適解よりも初期値に近い局所的最適解が存在する可能性があり、非線形最適化演算による解の更新を繰り返すうちに、解が大局的最適解に更新される前に、局所的最適解に更新される可能性が高いからである。従って、非線形最適化演算により局所的最適解が算出されることを防止するためには、初期値として、大局的最適解にできるだけ近い値を採用する必要がある。本実施形態は、初期楕円面VE0に基づいて、初期楕円面補正行列T0及び初期中心点cE0を算出し、これらの値を非線形最適化演算の初期値として適用する。複数の磁気データq1~qNの示す座標を近傍に有するように定められた初期楕円面VE0は、複数の磁気データq1~qNの示す座標との誤差を最小化する楕円面VEと近い形状を有する。従って、初期楕円面補正行列T0及び初期中心点cE0の示す座標は、大局的最適解として算出されるべき値である楕円面補正行列TE及び中心点cOGの示す座標に近い値である図7及び図12参照。このような初期値を用いた非線形最適化演算を行うことにより算出される最適楕円面補正行列TOP及び最適中心点cEOPの示す座標は、局所的最適解ではなく、大局的最適解厳密には、大局的最適解の近似値となる。このように、本実施形態の非線形最適化演算は、大局的最適解に近い値の適切な初期値が設定されるため、局所的最適解に陥ることなく、大局的最適解を最適解として算出することができる。ところで、本実施形態では、式69に示すように、非線形最適化演算において最適化の対象となる楕円面最適化関数fELT,cの変数として実対称行列である変数行列Tを用いるが、変数行列Tの代わりに、対称行列に限定されない実行列一般を表す変数行列TRを変数として用いて非線形最適化演算を行う方法も存在する非特許文献1参照。しかし、対称行列に限定されない変数行列TRは、任意のベクトルを変数行列TRの3つの固有ベクトル方向に各々伸縮させる座標変換に加え、任意のベクトルの向きを回転させる座標変換を表す場合がある。そのため、第1変数ベクトルqi-cを変数行列TRにより変換することで得られるベクトルTRqi-cは、第2磁気ベクトルTqi-cを任意角度回転させたベクトルとして算出される場合がある。すなわち、図14に示すように、変数行列TRを用いて算出される変換後磁気データROsiの示す座標は、変数行列Tを用いて算出される変換後磁気データsiの示す座標を、球面SEOP上で任意角度回転させた座標として算出される。この場合、変数行列TRを用いて算出される変換後磁気データROsiの示す座標に基づいて、地磁気Bgの向きを算出することは困難である。非特許文献1によれば、3次元磁気センサ60から見て方向が既知の参照磁界を用いることにより、変数行列TRが行う座標変換において生じた回転の角度を特定し、座標変換において回転が生じなかった場合の変換後磁気データsiの示す座標を算出している。この方法は、機器1が、参照磁界を測定する機会を有することが必要となる。これに対して、本実施形態では、変数行列Tを実対称行列に限定している。実対称行列は、互いに直交する3つの固有ベクトルと、3つの固有ベクトルに対応する3つの固有値とを有する。そして、実対称行列によりベクトルを変換する場合、変換後のベクトルは、変換前のベクトルを実対称行列の3つの固有ベクトル方向を向く3つのベクトルの和で表したとき、3つのベクトルの各々を、その方向を変えずに対応する固有値倍だけ伸縮させた3つのベクトルの和として算出される。すなわち、実対称行列は、任意のベクトルを、当該実対称行列が有する各固有ベクトル方向に各々伸縮させる座標変換を行う行列である。従って、実対称行列である変数行列Tを用いる非線形最適化演算は、最適楕円面補正行列TOPを、任意のベクトルを最適楕円面補正行列TOPの各固有ベクトル方向に各々伸縮させる座標変換を表す行列として算出するため、最適楕円面補正行列TOPは、回転を伴う座標変換は行わない。このような最適楕円面補正行列TOPにより、磁気データqiの示す座標を、変換後磁気データsiの示す座標へと変換することにより、地磁気Bgの正確な向きを得ることができる。また、変数行列TRは3行3列の行列で9個の独立した成分を有し、楕円面最適化関数fELTR,cの変数は12個となるのに対して、本実施形態に係る変数行列Tは対称行列であるため6個の独立した成分を有し、楕円面最適化関数fELT,cの変数は9個である。従って、本実施形態に係る非線形最適化演算は、変数行列TRを用いる非線形最適化演算に比べて変数の数が少なく、計算負荷が低減される。なお、楕円面最適補正値生成部400は、最適楕円面補正行列TOPが正定値行列であることを判定したうえで、最適楕円面補正行列TOP及び最適中心点cEOPを出力してもよい。最適楕円面補正行列TOPは、楕円面上の座標を、当該楕円面の主軸方向に伸縮させることにより、球面上の座標へと変換する行列であるため、最適楕円面補正行列TOPの有する3つの固有値は、いずれも正の値となるからである。楕円面最適補正値生成部400が出力した最適楕円面補正行列TOP及び最適中心点cEOPは、蓄積部100に格納される。[5.地磁気の算出]前述の通り、地磁気算出部600は、オフセット採用部610と、地磁気ベクトル計算部620とを備え、3次元磁気センサ60から出力される磁気データqiの示す座標に対して、楕円補正を施すことにより、地磁気SBgの向きを算出する図9参照。以下において、地磁気算出部600の動作を説明する。まず、オフセット採用部610は、楕円面補正行列TEとして最適楕円面補正行列TOPを採用するとともに、オフセットcOFFとして最適中心点cEOPの座標を示すベクトルを採用する。これにより、本実施形態に係る地磁気測定装置は、楕円補正を表す式1を式67に変形し、式67に基づいた楕円補正を行うことが可能となる。次に、地磁気ベクトル計算部620は、式67に基づいて楕円補正を行い、地磁気Bgの向きを算出する。具体的には、地磁気ベクトル計算部620は、オフセットcOFFである最適中心点cEOPの座標を起点とし、磁気データqiの示す座標を終点とする第1磁気ベクトルqi-cEOPを、楕円面補正行列TEである最適楕円面補正行列TOPにより変換し、第2磁気ベクトルsi-cEOPを算出する。このとき、第2磁気ベクトルsi-cEOPは、ずれ角ψを考慮しなければ、地磁気Bgと同じ向きを向く。地磁気ベクトル計算部620は、必要な場合にはずれ角ψを考慮して図7、及び段落0046参照、第2磁気ベクトルsi-cEOPから、地磁気Bgの向きを算出する。なお、本実施形態では、説明の便宜上、第2磁気ベクトルsi-cEOPを、最適中心点cEOPの示す座標を起点として配置しているが図7参照、第2磁気ベクトルsi-cEOPを、センサ座標系ΣSの原点を起点として配置しても構わない。この場合、球面SEOPは、センサ座標系ΣSの原点を中心とする半径1の球面を表し、変換後磁気データsiは原点を中心とする半径1の球面近傍に分布する。[6.第1実施形態の結論]以上に説明した第1実施形態では、複数の磁気データq1~qNの示す座標が、楕円面近傍に分布すると仮定して、複数の磁気データq1~qNの示す座標を近傍に有する楕円面を特定し、複数の磁気データq1~qNの示す座標を、当該楕円面と同一の中心点を有する球面近傍に座標変換する楕円補正を行った。これにより、3次元磁気センサ60を搭載する機器が軟磁性材料を備え、ソフトアイアン効果が生じている場合に、複数の磁気データq1~qNに基づいて、地磁気Bgの正確な方向を算出することができる。なお、3次元磁気センサ60を搭載する機器が軟磁性材料を備えず、ソフトアイアン効果が生じていない場合であっても、複数の磁気データq1~qNの示す座標が楕円面の近傍に分布する場合がある。例えば、3次元磁気センサ60が備える、X軸地磁気センサ61、Y軸地磁気センサ62、及びZ軸地磁気センサ63の3つセンサの感度にバラつきがある場合、センサ座標系ΣSにおいて本来球面近傍に分布するはずの複数の磁気データq1~qNの示す座標は、3つのセンサの感度に応じて当該球面をセンサ座標系ΣSの各軸方向に伸縮させた楕円面近傍に分布する。すなわち、センサ感度のバラつきがある場合、複数の磁気データq1~qNの示す座標は、センサ座標系ΣSの3つの軸の向きと同じ方向を向く3つの主軸を有する楕円面近傍に分布する。第1実施形態では、複数の磁気データq1~qNの示す座標が楕円面近傍に分布しているものと看做して楕円補正を行う。この楕円補正に用いる楕円面の3つの主軸は、直交関係を保ちつつ任意の方向を向くことができるため、第1実施形態に係る楕円補正は、楕円面の主軸とセンサ座標系ΣSの3軸の方向が一致する場合についても、適用することが可能である。従って、第1実施形態の楕円補正は、センサ感度のバラつきが原因で複数の磁気データq1~qNの示す座標が楕円面近傍に分布する場合、つまり、楕円面の主軸とセンサ座標系ΣSの3軸の方向が一致する場合であっても、地磁気Bgの正確な向きを算出することができる。また、第1実施形態では、初期楕円面VE0の算出に当たり、空間Ω上の第1評価軸ξ1、第2評価軸ξ2、及び第3評価軸ξ3の3つの異なる評価軸を用いて、3つの異なる楕円面を生成したうえで、センサ座標系ΣSにおいてこれらの3つの異なる楕円面の各々の中心点の相互間の距離が全て第1閾値Δc以下であるか否かを判定した。そして判定結果が肯定となる場合に限り、初期楕円面補正行列T0及び初期中心点cE0の座標を算出した。これにより、複数の磁気データq1~qNの示す座標の分布の形状から、楕円面の形状を特定することが困難な場合に、複数の磁気データq1~qNの示す座標の分布の形状とは異なる不正確な初期楕円面VE0に基づいた、不適切な初期楕円面補正行列T0及び初期中心点cE0の座標が生成されることを防止することが可能となり、不正確な値がオフセットとして採用されることを防止することが可能となった。また、第1実施形態では、楕円面VEに近い形状を有する初期楕円面VE0に基づいて定められる初期楕円面補正行列T0の各成分及び初期中心点cE0の座標を初期値とする非線形最適化演算により、最適楕円面補正行列TOP及び最適中心点cEOPの座標を算出した。非線形最適化は、初期値に、大局的最適解に近い値を採用することで、局所的最適解が最適解として算出される可能性を低減させ、大局的最適解が最適解として算出される可能性を高める。従って、本実施形態に係る非線形最適化演算は、局所的最適解が最適解として算出される可能性を低減し、不適切な最適楕円面補正行列TOPを用いた楕円補正による不正確な地磁気Bgの向きが算出される可能性を低減した。また、第1実施形態では、実対称行列である変数行列Tを変数とする楕円面最適化関数fELT,cを最小化することで、最適楕円面補正行列TOP及び最適中心点cEOPを算出した。これにより、最適楕円面補正行列TOPは、任意のベクトルを、最適楕円面補正行列TOPの3つの固有ベクトル方向に伸縮させる座標変換を行う行列として算出されるため、回転を伴う座標変換は行わない。すなわち、最適楕円面VEOP近傍に位置する磁気データqiの示す座標を楕円補正することで算出される変換後磁気データsiの示す座標は、最適中心点cEOPから見て、地磁気Bgと同じ方向を向く座標として求められるため、変換後磁気データsiに基づいて、正確な地磁気Bgの向きを算出することができる。<B.第2実施形態>以下、第2実施形態について説明する。[7.第2実施形態に係る地磁気測定装置の概要]第1実施形態では、3次元磁気センサ60が測定する対象の磁界を、地磁気Bg、内部磁界Bi、及び着磁磁界Bmに限定し、複数の磁気データq1~qNの示す座標が楕円面近傍に分布する場合を想定した。しかし、機器1の外部に磁界を生じる物体が存在する場合、物体の発する外部磁界Bxの影響により、複数の磁気データq1~qNの示す座標は、楕円面とは異なる形状に分布することがある。この場合には、複数の磁気データq1~qNの示す座標が楕円面近傍に分布すると仮定して算出された楕円面の中心点の示す座標は、3次元磁気センサ60のオフセットを示す座標とは一致しない可能性が高く、楕円面の中心点の示す座標をオフセットとして補正を行っても正確な地磁気Bgの向きを求めることはできない。また、機器1が軟磁性材料21を備えない場合、ソフトアイアン効果は生じないため、複数の磁気データq1~qNの示す座標は、図3に示したように球面近傍に分布し、楕円面近傍には分布しない。この場合、楕円補正を行わなくても、正確な地磁気Bgの向きを求めることができる。第2実施形態は、機器1の外部の物体の発する外部磁界Bxが存在する場合、及び、機器1が軟磁性材料21を備えず着磁磁界Bmが存在しない場合に対応した地磁気測定装置を実現することを目的とする。図15は、第2実施形態に係る地磁気測定装置のうち、CPU10が磁気データ処理プログラムを実行することにより実現される機能を表す、機能ブロック図である。第2実施形態に係る地磁気測定装置は、楕円面球面変換部500、分布判定部700、中心点算出部800、及び歪判定部900を備える点と、地磁気算出部600の代わりに地磁気算出部600aを備える点とを除き、第1実施形態に係る地磁気測定装置図9参照と同様に構成される。楕円面球面変換部500は、式67に基づいて、最適楕円面補正行列TOP、最適中心点cEOP、及び、複数の磁気データq1~qNから、複数の変換後磁気データs1~sNを算出する。具体的には、楕円面球面変換部500は、まず、式67に示ように、最適中心点cEOPの座標を起点として磁気データqiの示す座標を終点とする第1磁気ベクトルqi-cEOPを、最適楕円面補正行列TOPにより変換し、最適中心点cEOPを起点として変換後磁気データsiの示す座標を終点とする第2磁気ベクトルsi-cEOPを算出することで、変換後磁気データsiの示す座標を算出する。その後、楕円面球面変換部500は、算出した複数の変換後磁気データs1~sNを蓄積部100のバッファBU2に格納する。分布判定部700は、センサ座標系ΣSにおいて、複数の磁気データq1~qNの示す座標の分布が、3次元的な広がりを有しているか否かを判定し、判定結果を出力する。中心点算出部800は、センサ座標系ΣSにおいて、複数の磁気データq1~qNの示す座標を近傍に有する球面Sの中心点cSの示す座標を算出する。図3において説明したように、3次元磁気センサ60が測定する磁界が、地磁気Bg及び内部磁界Biである場合、複数の磁気データq1~qNの示す座標が、球面SGの近傍に分布するため、球面Sと球面SGとは一致すると看做すことができ、球面Sの中心点cSの示す座標が、オフセットcOFFを表す。第2実施形態に係る地磁気測定装置は、中心点算出部800を備えることで、着磁磁界Bmが存在せず、複数の磁気データq1~qNの示す座標が球面近傍に分布する場合に、簡易な方法により、3次元磁気センサ60のオフセットcOFFを算出することが可能となる。歪判定部900は、センサ座標系ΣSにおいて、複数の入力座標複数の磁気データq1~qNの示す座標、または複数の変換後磁気データs1~sNの示す座標がある立体近傍に分布すると仮定して、当該立体の形状と球面の形状との相違の程度を評価することで、立体の形状が球面と看做すことができる否かを判定し、判定結果を出力する。外部磁界Bxが存在し、複数の磁気データq1~qNの示す座標が球面または楕円面とは異なる歪んだ形状の立体近傍に分布する場合、3次元磁気センサ60のオフセットcOFFの正確な値を算出することは、困難である。第2実施形態に係る地磁気測定装置は、歪判定部900を備えることで、外部磁界Bxの影響が大きく、オフセットcOFFの算出が困難な場合に、不正確なオフセットcOFFが算出されることを防止し、不正確なオフセットを用いた補正処理により不正確な地磁気Bgが算出されることを防止する。地磁気算出部600aは、オフセット採用部610の代わりにオフセット採用部610aを備える点を除き、地磁気算出部600と同様に構成される。オフセット採用部610aは、オフセットcOFFとして、最適中心点cEOPの座標を示すベクトルまたは球面Sの中心点cSの座標を示すベクトルのいずれか一方を採用する。そして、オフセット採用部610aは、オフセットcOFFとして中心点cSの座標を示すベクトルを採用する場合、楕円面補正行列TEとして3行3列の単位行列Iを採用する。このとき、地磁気ベクトル計算部620は、オフセットcOFFである中心点cSの示す座標と、楕円面補正行列TEである単位行列Iとを用いて、式1に基づいた楕円補正を行い、地磁気Bgの向きを算出する。具体的には、楕円面補正行列TEは単位行列Iであるため、地磁気ベクトル計算部620は、ベクトルqi-cSを、地磁気SBgの向きを表すベクトルとして算出する。なお、単位行列Iを用いる楕円補正とは、式1からも明らかなとおり、単に、磁気データqiの示す座標からオフセットcOFFを減算する演算に過ぎず、実質的には楕円補正を行わないことを意味する。従って、楕円面補正行列TEとして単位行列Iが採用された場合、地磁気ベクトル計算部620は、式1に基づいた演算を行わず、単純に、磁気データqiの示す座標からオフセットcOFFとして採用された中心点cSの座標を減算する処理を行ってもよい。一方、オフセット採用部610aは、オフセットcOFFとして最適中心点cEOPの座標を示すベクトルを採用する場合、楕円面補正行列TEとして最適楕円面補正行列TOPを採用する。このとき、地磁気ベクトル計算部620は、オフセットcOFFである最適中心点cEOPの示す座標と、楕円面補正行列TEである最適楕円面補正行列TOPとを用いて、式1に基づいた楕円補正を行い、地磁気Bgの向きを算出する。具体的には、地磁気ベクトル計算部620は、式1を変形した式67により、第2磁気ベクトルsi-cEOPを、地磁気SBgの向きを表すベクトルとして算出する。以下では、第2実施形態に係る地磁気測定の処理の詳細を説明する前提として、外部磁界Bxの性質について説明する。図16は、3次元磁気センサ60の検出対象である、地磁気Bg、内部磁界Bi、着磁磁界Bm、及び外部磁界Bxを、地上座標系ΣGにおいて表した概念図である。ここで、図16に示す位置PSは、地上座標系ΣGにおけるセンサ座標系ΣSの原点の位置すなわち、地上座標系ΣGにおける、3次元磁気センサ60の位置を表す。第2実施形態に係る地磁気測定装置は、機器1の他に、図16に示す機器1aに適用することができる。ここで、機器1aは、部品2の代わりに、軟磁性材料21を含まない部品2aを備える点を除き、機器1と同様の構成である。すなわち、機器1aは、機器1と異なり、着磁磁界Bmを発しない。図16に示すように、外部磁界Bxは、機器1または機器1aの外部に存在する物体3から発せられる磁界であり、物体3との相対的位置関係によって向き及び大きさが変化する不均一な磁界である。地上座標系ΣGにおける3次元磁気センサ60の位置Psを変化させた場合、3次元磁気センサ60が測定する外部磁界Bxの向き及び大きさが変化する。従って、外部磁界Bxは、地上座標系ΣGにおいて、位置Psに依存して方向及び大きさの双方を変化させるベクトルGBxPsとして表現される。また、地上座標系ΣGにおける3次元磁気センサ60の姿勢μを変化させた場合、3次元磁気センサ60が測定する外部磁界Bxの向きが変化する。図17は、3次元磁気センサ60の位置PsをPS1~PSNと変化させると共に、姿勢μをμ1~μNと変化させて磁界を測定したときの、3次元磁気センサ60が出力する複数の磁気データq1~qNを、センサ座標系ΣSにおいてプロットした図である。なお、図17では、簡単のために、着磁磁界Bmは存在せず、内部磁界Bi、地磁気Bg、及び外部磁界Bxのみが存在する場合を想定する。外部磁界Bxは、センサ座標系ΣSにおいて、3次元磁気センサ60の位置Psに依存して方向及び大きさの双方を変化させ、3次元磁気センサ60の姿勢μに依存して方向を変化させるベクトルSBxμ、Psとして表現される。3次元磁気センサ60が、内部磁界Bi、地磁気Bg、及び外部磁界Bxを測定する場合、複数の磁気データq1~qNの各々の表す座標は、内部磁界を表すベクトルSBi、地磁気を表すベクトルSBgμ、及び外部磁界を表すベクトルSBxμ、Psの和を表すベクトルにより示される。従って、複数の磁気データq1~qNの示す座標は、中心点cOGを起点とする地磁気を表すベクトルSBgμの終点を表す球面SGと、中心点cOGを起点とする外部磁界を表すベクトルSBxμ、Psの終点を表す曲面SXとを、中心点cOGを起点として重ね合わせた立体SDの表面近傍に分布する。外部磁界Bxを表す曲面SXが、球面とは異なる歪んだ形状を有する場合には、立体SDも球面とは異なる歪んだ形状を有する。立体SDが、球面とは異なる歪んだ形状を有する場合、複数の磁気データq1~qNの示す座標に基づいて、地磁気Bgを表す球面SGの中心点cOGの座標を算出することは困難である。仮に、複数の磁気データq1~qNの示す座標を近傍に有するような球面Sを定め、球面Sの中心点cSを算出しても、球面Sの中心点cSと球面SGの中心点cOGとは異なる座標となる可能性が高いからである。図24参照。このように、立体SDが球面とは異なる歪んだ形状を有し、球面SGの中心点cOGの示す座標の算出が困難である場合には、複数の磁気データq1~qNに基づいたオフセットcOFFの算出を防止する必要がある。但し、不均一な外部磁界Bxの影響が小さく、立体SDの形状がほぼ球面と看做せる場合には、複数の磁気データq1~qNの示す座標に基づいて、球面SGの中心点cOGの示す座標を算出することができる。例えば、図18Aに示すように、外部磁界Bxが微弱である場合、地磁気Bgを表す球面SGと、外部磁界Bxを表す曲面SXとの重ね合わせである立体SDは、球面SGとほぼ同一の図形となる。従って、複数の磁気データq1~qNの示す座標は、球面SGの近傍に分布すると看做すことができるため、複数の磁気データq1~qNから、球面SGの中心点cOGの示す座標を算出することができる。また、図18Bに示すように、不均一な外部磁界Bxが大きい場合であっても、立体SDの形状がほぼ球面と看做せる場合がある。例えば、不均一な外部磁界Bxが存在する場合であっても、N個の磁気データq1~qNを取得する際に、機器1または機器1aの利用者が機器1または機器1aを手で握って3次元磁気センサ60の位置PSが変化するように振るのではなく、3次元磁気センサ60の位置PSを固定して姿勢μのみを変化させる場合には、外部磁界Bxは、センサ座標系ΣSにおいて、姿勢μに基づいてその方向のみを変化させる一定の大きさのベクトルSBxμとして表現される。この場合、外部磁界Bxを表す曲面SXの形状は、中心点cOGを中心とする球面となるため、中心点cOGを中心とする球面SGと、中心点cOGを中心とする球面の形状の曲面SXとを、中心点cOGを中心に重ね合わせた立体SDの形状は、中心点cOGを中心とする球面となる。従って、複数の磁気データq1~qNの示す座標に基づいて、立体SDの表す球面の中心点の座標を算出することにより、球面SGの中心点cOGの示す座標を算出することができる。本実施形態では、複数の磁気データq1~qNの示す座標に基づいて、外部磁界Bxによる影響の大きさ、すなわち、立体SDの形状と球面の形状との相違の程度を評価する。これにより、複数の磁気データq1~qNの示す座標に基づいて、オフセットcOFFが算出可能であるかを判断し、外部磁界Bxの影響を受けた不正確なオフセットが算出されることを防止することが可能となる。なお、詳細は後述するが、本実施形態に係る地磁気測定装置は、楕円面補正部200、楕円面球面変換部500、及び歪判定部900により、複数の磁気データq1~qNの示す座標の分布の形状と、楕円面の形状との相違の程度について評価することもできる。複数の磁気データq1~qNの示す座標を、楕円補正により複数の変換後磁気データs1~sNの表す座標に変換した場合、複数の変換後磁気データs1~sNの示す座標を近傍に有する立体SDEの形状が球面と看做すことができれば、複数の磁気データq1~qNの示す座標を近傍に有する立体SDの形は、楕円面と看做すことができるからである。すなわち、図15に示すように、楕円面補正部200、楕円面球面変換部500、及び歪判定部900は、複数の磁気データq1~qNの示す座標の分布の形状が、球面と看做すことのできる形状、楕円面と看做すことのできる形状、または、球面と看做すことも楕円面と看做すこともできない歪んだ形状、のいずれに該当するかを判断する、歪形状判断部4として機能する。そして、歪形状判断部4が、複数の磁気データq1~qNの示す座標の分布の形状を、球面と看做すことのできる形状、または、楕円面と看做すことのできる形状であると判断した場合、地磁気測定装置は、オフセットcOFFを算出するが、球面とも楕円面とも異なる歪んだ形状であると判断した場合、オフセットcOFFを算出しない。このように、本実施形態に係る地磁気測定装置は、歪形状判断部4を備えることにより、外部磁界Bxの影響による不正確なオフセットが算出されることを防止することが可能になるとともに、外部磁界Bxの影響が無視できる程度に小さい場合には、ソフトアイアン効果が生じている場合及び生じていない場合の双方において正確なオフセットcOFFの算出が可能になる。以下、本実施形態における、オフセットcOFFの候補となる座標の算出方法と、当該座標のオフセットcOFFへの採否の判断方法について、具体的に説明する。[8.第2実施形態に係る地磁気測定装置のオフセット導出処理フロー]図19は、第2実施形態に係る地磁気測定装置が、オフセットを導出する動作を説明するためのフローチャートである。このフローチャートは、CPU10が、本実施形態に係る磁気データ処理プログラムを実行することにより実施される。ステップS1において、地磁気測定装置は、初期化処理を行う。初期化処理は、蓄積部100に備えられるバッファBU1に格納された複数の磁気データq1~qNと、蓄積部100に備えられるバッファBU2に格納された各種データ複数の変換後磁気データs1~sN等とを、破棄する処理である。なお、本実施形態の地磁気測定装置は、初期化処理において、バッファBU1に格納された複数の磁気データq1~qNの全部を廃棄するが、古い方から一定割合のデータのみを廃棄してもよい。ステップS2において、地磁気測定装置は、磁気データ取得処理を行う。磁気データ取得処理は、3次元磁気センサ60から順次出力される複数の磁気データq1~qNを、蓄積部100に備えられるバッファBU1に格納する処理であるNは精度よくオフセットを導出するために必要な磁気データの規定測定回数を表す9以上の自然数。ステップS3において、地磁気測定装置は、磁気データ分布判定処理を行う。磁気データ分布判定処理は、分布判定部700が行う処理である。磁気データ分布判定処理において、分布判定部700は、センサ座標系ΣSにおいて、複数の磁気データq1~qNの示す座標の分布が、3次元的な広がりを有しているか否かを判定し、判定結果を出力する。判定結果が肯定である場合、地磁気測定装置は、処理をステップS4に進める。一方、判定結果が否定である場合、すなわち、複数の磁気データq1~qNの示す座標の分布が2次元的または1次元的である場合、地磁気測定装置は、処理をステップS1に戻す。ステップS4において、地磁気測定装置は、中心点算出処理を行う。中心点算出処理は中心点算出部800が行う処理である。中心点算出処理において、中心点算出部800は、センサ座標系ΣSにおいて、複数の磁気データq1~qNの示す座標を近傍に有する球面Sの中心点cSの示す座標を算出し、これを出力する。ステップS5において、地磁気測定装置は、歪判定処理を行う。ステップS5における歪判定処理は、歪判定部900が行う処理である。ステップS5の歪判定処理において、歪判定部900は、複数の入力座標として複数の磁気データq1~qNの示す座標を適用したうえで、複数の入力座標がある立体SDの近傍に分布すると仮定して、立体SDの形状と球面の形状との相違の程度を評価することで、立体SDの形状が球面と看做すことができるか否かを判定し、判定結果を出力する。判定結果が肯定である場合、地磁気測定装置は、処理をステップS10に進める。一方、判定結果が否定である場合、地磁気測定装置は、処理をステップS6に進める。ステップS6において、地磁気測定装置は、初期楕円面生成処理を行う。初期楕円面生成処理は、3節において説明した楕円面初期補正値生成部300が行う処理である。前述の通り、初期楕円面生成処理において、楕円面初期補正値生成部300は、複数の磁気データq1~qNの示す座標を近傍に有する初期楕円面VE0の中心点である初期中心点cE0の座標と、初期楕円面VE0上の座標を球面SE0上の座標へと変換する初期楕円面補正行列T0とを算出する。また、楕円面初期補正値生成部300は、複数の磁気データq1~qNに基づいて算出される、第1楕円面係数行列Dxx、第2楕円面係数行列Dyy、第3楕円面係数行列Dzz、中心点cxxの座標、中心点cyyの座標、及び中心点czzの座標が、第1の条件及び第2の条件の双方を満たすか否かの判定を行う。当該判定結果が肯定である場合、地磁気測定装置は、処理をステップS7に進める。一方、判定結果が否定である場合、地磁気測定装置は、処理をステップS1に戻す。なお、前述のとおり、楕円面初期補正値生成部300は、第1の条件についての判定を行わず、第2の条件についてのみ判定を行うものであってもよい。ステップS7において、地磁気測定装置は最適楕円面生成処理を行う。最適楕円面生成処理は、4節において説明した楕円面最適補正値生成部400の行う処理である。前述の通り、最適楕円面生成処理において、楕円面最適補正値生成部400は、初期楕円面補正行列T0及び初期中心点cE0に基づいて最適楕円面補正行列TOP及び最適中心点cEOPの座標を算出する。ステップS8において、地磁気測定装置は楕円面球面変換処理を行う。楕円面球面変換処理は、楕円面球面変換部500が行う処理である。楕円面球面変換処理において、楕円面球面変換部500は、最適楕円面補正行列TOP及び最適中心点cEOPに基づいて、最適楕円面VEOP近傍に存在する複数の磁気データq1~qNの示す座標を、複数の変換後磁気データs1~sNにより表される球面SEOP近傍の座標へと変換する。その後、楕円面球面変換部500は、複数の変換後磁気データs1~sNを、蓄積部100に備えられたバッファBU2に格納する。バッファBU2は、RAM20により構成される。ステップS9において、地磁気測定装置は、歪判定処理を行う。ステップS9における歪判定処理は、ステップS5における歪判定処理と同様に、歪判定部900が行う処理である。ステップS9の歪判定処理において、歪判定部900は、複数の入力座標として複数の変換後磁気データs1~sNの示す座標を適用したうえで、複数の入力座標が立体SDEの近傍に分布すると仮定して、立体SDEの形状と球面の形状との相違の程度を評価することで、立体SDEの形状が球面と看做すことができるか否かを判定し、判定結果を出力する。判定結果が肯定である場合、地磁気測定装置は、処理をステップS10に進める。一方、判定結果が否定である場合、地磁気測定装置は、処理をステップS1に戻す。なお、以下では、ステップS5において行われる歪判定処理と、ステップS9において行われる歪判定処理とを区別する場合、前者を第1の歪判定処理と称し、後者を第2の歪判定処理と称する。また、第1の歪判定処理において形状が評価される立体SDと区別するため、第2の歪判定処理において形状が評価される立体を立体SDEと称する。第1の歪判定処理と第2の歪判定処理とは、複数の入力座標の値が異なる以外は、同一の処理である。ステップS10において、地磁気測定装置は、オフセット採用処理を行う。オフセット採用処理は、オフセット採用部610aが行う処理である。オフセット採用処理において、オフセット採用部610aは、中心点cSまたは中心点cEOPの示す座標をオフセットcOFFとして採用するとともに、単位行列Iまたは最適楕円面補正行列TOPを楕円面補正行列TEとして採用する。具体的には、ステップS5における判定結果が肯定である場合、オフセット採用部610aは、ステップS4において中心点算出部800が算出した球面Sの中心点cSの座標を表すベクトルをオフセットcOFFとして採用するとともに、単位行列Iを楕円面補正行列TEとして採用する。一方、ステップS5における判定結果が否定であり、且つ、ステップS9における判定結果が肯定である場合、オフセット採用部610aは、ステップS7で楕円面最適補正値生成部400が算出した最適中心点cEOPの座標を表すベクトルをオフセットcOFFとして採用するとともに、楕円面最適補正値生成部400が算出した最適楕円面補正行列TOPを楕円面補正行列TEとして採用する。そして、オフセット採用部610aは、オフセットcOFFと楕円面補正行列TEとを出力する。また、ステップS9における判定結果が否定である場合、オフセット採用部610aは、オフセットcOFF及び楕円面補正行列TEを採用しない。なお、5節において説明したとおり、地磁気ベクトル計算部620は、オフセットcOFFと楕円面補正行列TEとを用いて、3次元磁気センサ60から出力される磁気データqiの示す座標に対して楕円補正を施すことで、地磁気Bgの向きを算出する。地磁気ベクトル計算部620が楕円補正に用いるオフセットcOFF及び楕円面補正行列TEは、オフセット採用部610aから出力されるオフセットcOFF及び楕円面補正行列TEにより更新される。本実施形態では、ステップS9において、判定結果が否定である場合、地磁気測定装置は、処理をステップS1に戻すが、その際に、何らかのメッセージを表示部50に出力した上で処理を一旦停止させ、ユーザからの指示を待ってステップS1の処理を再開させてもよい。N個の磁気データq1~qNを取得する際に、利用者が機器1を手で握って回転させるのではなく、機器1の位置を固定したままその姿勢のみを変化させるようにすると、外部磁界Bxの影響を低く抑えることができる図18B参照。そのため、ステップS9における判定結果が否定である場合には、機器1の位置を固定したまま回転させることを、利用者に対して指示してもよい。利用者に対する指示は、機器1の表示部50に画像や動画等を表示したり、音声等を用いたりすることで、行ってもよい。また、本実施形態では、ステップS6またはS9において、判定結果が否定である場合に、処理をステップS1に戻さずに、フローチャートに示された処理を終了してもよい。このように、本実施形態の地磁気測定装置は、ステップS5において、複数の磁気データq1~qNの示す座標が、球面Sの近傍に分布するか、球面とは異なる歪んだ形状の立体SDの近傍に分布するかを判定する。そして、複数の磁気データq1~qNの示す座標が、球面とは異なる歪んだ形状の立体SDの近傍に分布すると判定された場合、ステップS6~S9を実行することにより、複数の磁気データq1~qNの示す座標が、最適楕円面VEOPの近傍に存在するか否かの判定を行う。すなわち、本実施形態の地磁気測定装置は、複数の磁気データq1~qNの示す座標の分布形状が、球面、楕円面、及び球面とも楕円面とも異なる歪んだ形状の立体のうち、いずれの形状と看做すことが適当であるかについて、判別することが可能となる。これにより、本実施形態に係る地磁気測定装置は、複数の磁気データq1~qNの示す座標が、球面または楕円面の近傍に分布する場合には、これらの中心点をオフセットとして採用することにより、正確な地磁気の向きを算出することができる。一方、複数の磁気データq1~qNの示す座標が、球面とも楕円面とも異なる歪んだ形状の立体の近傍に分布する場合、本実施形態に係る地磁気測定装置は、オフセットの算出を防止するため、不正確な地磁気の向きの算出を防止することができる。以下において、磁気データ分布判定処理、中心点算出処理、及び歪判定処理の詳細を説明する。なお、理解を容易なものとするために、中心点算出処理を説明した後に、磁気データ分布判定処理の説明を行う。[9.中心点算出処理]図20を参照しつつ、ステップS4において、中心点算出部800が行う中心点算出処理について説明する。中心点算出処理は、3次元磁気センサ60が出力するN個の磁気データq1~qNの示す座標が、半径rSの球面Sの近傍に分布すると仮定して、球面Sの中心点cSの座標を算出する処理である。球面Sは、センサ座標系ΣS上で、複数の磁気データq1~qNの示す座標を近傍に有するように定められる球面の中心点の座標を求めるために、計算の便宜上導入される球面であり、地磁気Bgを表す球面SGとは異なる球面である。なお、以下において記載されるベクトル及び座標は、特に断りがない場合には、センサ座標系ΣSにおいて表現されたものであるとする。複数の磁気データq1~qNの示す座標を近傍に有する球面Sの中心点cSの座標の算出は、公知の方法を適宜適用すればよい。例えば、以下の方法により算出してもよい。磁気データqiの示す座標を式11で表し、中心点cSの座標を以下の式72で表した場合、複数の磁気データq1~qNの示す座標が、半径rSの球面S上に存在することは、以下の式71で表される。図20に示すように、磁気データqiの示す座標を、複数の磁気データq1~qNの重心qCの示す座標を起点とするベクトルqi-qCで表現した場合、式71に対して複数の磁気データq1~qNの示す座標の各々を代入した式と、式71とに基づいて、以下の式73を得ることができる。以下では、式73に示す方程式を、球面方程式と呼ぶ。ここで、重心qCは、以下の式74及び式75で定義される3次元のベクトルである。また、行列Xは、式76に示すN行3列の行列であり、ベクトルjは、式77に示すN次元のベクトルであり、値RAVEは、式78に示す値である。式73に示した球面方程式は、複数の磁気データq1~qNの示す座標の全てが、中心点cSを中心とする球面S上に完全に一致する場合には解を有する。しかし、3次元磁気センサ60の測定誤差等を考慮すると、複数の磁気データq1~qNの全てが、球面Sと完全に一致することはないため、球面方程式は解を持たない。そこで、統計的な手法により、球面方程式の尤もらしい解を得るために、式79で表す誤差を吸収するベクトルである第1球面誤差ベクトルδSを導入する。ここで、式79に現れる変数ベクトルcは、式5に示した3次元のベクトルであるが、本節では、中心点cSの座標を表すための変数として用いる。第1球面誤差ベクトルδSのノルムを最小にするベクトルc、換言すれば、Tを最小化するようなベクトルcの表す座標が、球面Sの中心点cSの座標として尤もらしいものであるといえる。ここで、以下の式80で示される中心点算出関数fScを定義すると、中心点算出関数fScを最小化するベクトルcの表す座標が、球面Sの中心点cSの座標として尤もらしい値となる。中心点cSの座標は、式82に示す3行3列の共分散行列Aが正則である場合には、式81により算出される。前述の通り、3次元磁気センサ60が、内部磁界Bi及び地磁気Bgのみを検出する場合には、球面Sと、地磁気Bgを表す球面SGとは、ほぼ同一の球面となり、球面Sの中心点cSと、球面SGの中心点cOGとは、ほぼ同一の座標となる。従って、3次元磁気センサ60が内部磁界Bi及び地磁気Bgのみを検出する場合、式81に示される中心点cSの座標を示すベクトルを、磁気センサのオフセットcOFFとして採用することができる。[10.磁気データ分布判定処理]以下では、ステップS3において、分布判定部700が行う磁気データ分布判定処理について説明する。前述した、中心点算出処理において、球面Sの中心点cSを算出するためには、複数の磁気データq1~qNの示す座標が、センサ座標系ΣSにおいて3次元的な広がりを有して分布していることを必要とする。しかし、機器13次元磁気センサ60の姿勢μは、機器1の利用者が手で機器1を動かすことにより変化するため、機器1の動かし方が不十分な場合には、機器1の姿勢変化は3次元的とはならず2次元的なものとなることがある。この場合、センサ座標系ΣSにおける複数の磁気データq1~qNの示す座標は、3次元的な広がりを有さず2次元的に分布する。例えば、図21に示すように、複数の磁気データq1~qNの示す座標が、センサ座標系ΣSの平面π上の円πC近傍に2次元的に分布する場合、球面Sは、円πCを切断面に有するような球面であるということしか特定できない。円πCを切断面に有するような球面は、円πCの中心点πCOを通り平面πに直交する直線πL上の中心点cπ1を中心とする球面Sπ1であるかもしれないし、直線πL上の中心点cπ2を中心とする球面Sπ2であるかもしれない。つまり、球面Sの中心点cSは、直線πL上に位置することまでは特定可能であるが、具体的に直線πL上のどの位置に存在するかについては特定不可能である。このように、複数の磁気データq1~qNの示す座標が2次元的に分布する場合、複数の磁気データq1~qNに基づいて、正確な中心点cSを算出することはできない。複数の磁気データq1~qNに基づいて球面Sの中心点cSを算出するためには、図22に示すように、センサ座標系ΣSにおいて複数の磁気データq1~qNの示す座標が3次元的に広がりを持って分布していることが必要となる。磁気データ分布判定処理において、分布判定部700は、複数の磁気データq1~qNの示す座標が3次元的に分布しているか否かを判定する。複数の磁気データq1~qNの示す座標が3次元的に分布しているか否かの判定は、公知の方法を適宜適用すればよいが、例えば、式82に示される共分散行列Aを用いて判定してもよい。以下において、共分散行列Aの性質を説明する。共分散行列Aの固有値を、大きい順番に最大固有値λ1、中間固有値λ2、最小固有値λ3とし、それぞれの固有値に対応する大きさ1に正規化された固有ベクトルをu1、u2、u3とする。また、磁気データqiを、前述した重心qCを原点とする重心座標系ΣCにおいて表したベクトルをCqiと表す。このとき、固有値λjj=1、2、3は、固有ベクトルuj方向の分散σ2jに等しい。図22に示すように、各固有ベクトルu1~u3を重心座標系ΣCの原点qCを起点となるように配置する。このとき、例えばj=1の場合について検討する。固有値λ1は、ベクトルCqiを、固有ベクトルu1へ射影した長さLi1の二乗Li12を、N個の磁気データCqii=1、2、...Nについて平均した値に等しい。すなわち、固有値λjは、N個の磁気データCqiが、重心qCから固有ベクトルuj方向にどの程度離れているか、つまり、複数の磁気データq1~qNの示す座標の分布が固有ベクトルujの方向にどの程度の広がり有するかを表す。最小固有値λ3に対応する固有ベクトルu3の方向が、複数の磁気データq1~qNの示す座標の分布の広がりが最も小さい方向であり、最小固有値λ3が、複数の磁気データq1~qNの示す座標の分布の広がりが最も小さい方向における広がりの程度を示す指標である。従って、複数の磁気データq1~qNの示す座標が3次元的に分布しているといえるためには、最小固有値λ3の値が、一定の閾値分散許容値λ0以上であればよい。磁気データ分布判定処理において、分布判定部700は、共分散行列Aの最小固有値λ3が閾値λ0以上であれば、複数の磁気データq1~qNの示す座標の分布が十分に3次元的であると判断して、処理を前述したステップS4の中心点算出処理へと進める。一方、分布判定部700は、最小固有値λ3が閾値λ0未満である場合には、複数の磁気データq1~qNの示す座標が3次元的な広がりを有さないと判断し、処理をステップS1の初期化処理に戻す。[11.歪判定処理]歪判定部900は、ステップS5において第1の歪判定処理を行うとともに、ステップS9において第2の歪判定処理を行う。第2の歪判定処理は、複数の入力座標として、複数の磁気データq1~qNの示す座標の代わりに、複数の変換後磁気データs1~sNの示す座標を用いて歪判定処理を行う点を除いて、第1の歪判定処理と同様の処理である。以下11.1節において、第1の歪判定処理について述べた上で、11.2節において、第2の歪判定処理について述べる。[11.1.第1の歪判定処理]歪判定処理は、複数の入力座標、すなわち、複数の磁気データq1~qNの示す座標が、球面とは異なる歪んだ形状の立体SDの表面近傍に分布すると仮定する。立体SDは、図23に示すように、球面第2球面S2と、歪誤差ベクトルkEとを、足し合わせることで得られる図形であり、以下の式83により表される。以下では、式83に示す方程式を、立体方程式と呼ぶ。ここで、球面S2は、中心点第2球面の中心点cS2を中心とする球面であり、立体方程式のうち、歪誤差ベクトルkE以外の成分Xc-qC-jとして表される。歪誤差ベクトルkEは、以下の式84により表されるN次元ベクトルである。但し、歪評価行列Eは、以下の式85に表される3行3列の対称行列であり、基準点wKEは、以下の式86により表される3次元ベクトルである。また、式83の右辺に現れる0NはN次元の零ベクトルである。式83の左辺に現れる変数ベクトルcは、式5で示された3次元のベクトルであるが、本節では、球面S2の中心点cS2を表すための変数として用いる。歪判定処理は、立体方程式のうち、歪を表す成分kEの大きさを評価することにより、立体SDの形状と、球面S2の形状とが、どの程度相違しているかについて評価する。具体的には、立体方程式における歪誤差ベクトルkEの影響の大きさを、後述する式93及び式94で示される歪評価値gDEに基づいて評価する。N次元の歪誤差ベクトルkEを構成するN個の要素のうち、第i行目の要素keは、以下の式87で表される関数kevに対して、基準点wKEの示す座標を起点として磁気データqiの示す座標を表すベクトルqi-wKEを代入することにより与えられる。関数kevは、式85に示す歪評価行列Eを係数行列とし、式88に示すベクトルvの3つの要素を変数とする、二次形式で表現される関数である。つまり、関数kevは、ベクトルvと、ベクトルvを歪評価行列Eにより変換したベクトルEvとの内積を示す。なお、第1の歪判定処理では、基準点wKEとして、以下の式89に示すように、球面Sの中心点cSを採用する。3次元磁気センサ60の測定誤差等を考慮すると、複数の磁気データq1~qNの示す座標の全てが、立体SDと完全に一致する位置に存在することはないため、式83に示した立体方程式は解を持たない。そこで、統計的な手法により、立体方程式の解として尤もらしい値を得るために、式90で表される誤差を吸収するベクトルである立体誤差ベクトルδSDを導入する。立体誤差ベクトルδSDは、以下の式91に示す第2球面誤差ベクトルδS2と、歪誤差ベクトルkEとを加算したものである。第2球面誤差ベクトルδS2は、立体方程式のうち、球面S2を表す成分Xc-qC-jに対応する成分である。立体誤差ベクトルδSDは、複数の磁気データq1~qNの示す座標と、立体SDの表面との誤差を表す、N次元のベクトルである。立体誤差ベクトルδSDのノルムを最小にする変数ベクトルc及び歪評価行列E、換言すれば、以下の式92で示す歪評価関数fSDE、cを最小化する変数ベクトルc及び歪評価行列Eによって、複数の磁気データq1~qNの示す座標を表面近傍に有する立体SDが表現される。以下で、式90に示した立体誤差ベクトルδSDの性質を、式79に示した第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の中心点cSを求める処理である。一方、立体誤差ベクトルδSDは、第2球面誤差ベクトルδS2と歪誤差ベクトルkEとの和によって表されるベクトルであり、磁気データq1~qNの示す座標と立体SDとの誤差を吸収するベクトルである。第2球面誤差ベクトルδS2は、第1球面誤差ベクトルδSと同様に、磁気データq1~qNの示す座標と、球面S2との誤差を、ホワイトノイズとして表現するベクトルである。一方、歪誤差ベクトルkEは、式87で示した3変数二次形式の関数kevを各要素とするベクトルである。3変数の二次形式は、変数が2次の項から構成される関数であり、3次元空間上の様々な曲面、例えば、直線、平面、柱面、球面、楕円面、錐面、1葉双曲面、2葉双曲面、及び各種放物面等を描くことができる。従って、歪誤差ベクトルkEは、各磁気データq1~qNが示す座標と球面S2とのN個の誤差の各々を、互いに独立したものではなく、N個の誤差の全てが同一の関数kevにより表される3次元空間上の曲面上に存在するという制約を持った値として表現する。このように、立体誤差ベクトルδSDは、各磁気データq1~qNが示す座標と球面S2とのN個の誤差を、ホワイトノイズである第2球面誤差ベクトルδS2と、球面S2からの歪を表す曲面を表す歪誤差ベクトルkEとに分離して表現する。立体方程式における歪誤差ベクトルkEの影響が、無視できる程度に小さい場合、立体SDと球面S2とは、等しい図形と看做すことができ、式92に定める歪評価関数fSDE、cと、式80に定める中心点算出関数fScとは、等しい関数と看做すことができる。このとき、歪評価関数fSDE、cを最小化して求められる立体SDと、中心点算出関数fScを最小化して求められる球面Sとは、等しいものと看做すことができるため、立体SDの表面近傍に分布する複数の磁気データq1~qNが示す座標は、球面Sの近傍にも分布すると看做すことができる。前述の通り、立体SDの形状が球面と看做せる場合、立体SDが表す球面の中心点の示す座標と、球面SGの中心点cOGとは、一致すると看做すことができる。よって、球面Sの中心点cSの示す座標と、球面SGの中心点cOGの示す座標は、等しいと看做すことができる。このように、立体方程式における歪誤差ベクトルkEの影響が小さい場合、中心点算出手段において算出した球面Sの中心点cSの示す座標は、球面SGの中心点cOGの示す座標と等しいと看做されるため、中心点cSの示す座標をオフセットcOFFとして採用することができる。逆に、立体方程式における歪誤差ベクトルkEの影響が大きい場合、図23に示すように、複数の磁気データq1~qNが示す座標と球面S2との誤差は、ホワイトノイズである第2球面誤差ベクトルδS2と、球面S2からの歪を表す歪誤差ベクトルkEとによって吸収される。この場合、立体SDの形状は、球面とは異なる形状を有する。また、立体方程式における歪誤差ベクトルkEの影響が大きい場合、歪評価関数fSDE、cと、中心点算出関数fScとは、異なる関数となる。この場合、図24に示すように、歪評価関数fSDE、cを最小化して求められる立体SDと、中心点算出関数fScを最小化して求められる球面Sとは、異なる図形となり、立体SDの表面近傍に分布する複数の磁気データq1~qNの示す座標は、球面Sの近傍に分布すると看做すことはできない。中心点算出処理は、複数の磁気データq1~qNの示す座標が球面Sの近傍に存在していることを前提として、球面SGの中心点cOGと等しいと看做すことのできる中心点cSの示す座標を算出する処理である。従って、複数の磁気データq1~qNの示す座標が球面Sの近傍に存在しない場合、中心点cSと、中心点cOGとは、一致しない。この場合、中心点cSの座標を示すベクトルを、オフセットcOFFとして採用することはできない。このように、立体方程式における歪誤差ベクトルkEの影響の大きさを評価することにより、球面Sの中心点cSを、オフセットcOFFとして採用することの可否を判断することが可能である。以下では、立体方程式おける歪誤差ベクトルkEの影響の大きさを評価する方法について述べる。ここで、以下の式93及び94に示す歪評価値gDEを、立体方程式における歪誤差ベクトルkEの影響の大きさを評価する評価値として定義する。歪評価値gDEは、歪評価行列Eの有する3つの固有値のうち、絶対値が最大となる最大固有値λE1の絶対値つまり、歪評価行列Eのノルムである。歪評価値gDEが、歪許容値δ0以下の小さな値であれば、立体SDと球面S2とは等しい図形であると看做すことができ、立体SDの表面近傍に分布する複数の磁気データq1~qNの示す座標は、球面Sの近傍に分布すると看做すこともできる。この場合、中心点算出処理により求められた球面Sの中心点cSの座標を表すベクトルを磁気センサのオフセットcOFFとして採用することができる。ところで、前述の通り、歪誤差ベクトルkEの各要素は、基準点wKEの示す座標から見た磁気データqiの示す座標を表すベクトルqi-wKEと、当該ベクトルqi-wKEを歪評価行列Eにより変換したベクトルEqi-wKEとの内積である。つまり、歪誤差ベクトルkEを構成する複数の要素の絶対値は、当該要素に対応する磁気データqiの示す座標を基準点wKEの示す座標から表したベクトルqi-wKEと、歪評価行列Eが有する3つの固有値のうち絶対値の大きさが最大となる最大固有値λE1に対応する固有ベクトルuE1とが平行となる場合には、大きな値となる。従って、球面S2との誤差が大きな座標を示す磁気データqiが多く存在する領域を、基準点wKEの示す座標から表した方向と、歪評価行列Eの最大固有値λE1に対応する固有ベクトルuE1の方向とが等しくなるように歪評価行列Eの各成分を定める場合、歪誤差ベクトルkEは、複数の磁気データq1~qNの示す座標と、球面S2との誤差の大きさを正確に表現する。歪評価関数fSDE、cを最小化する歪評価行列Eは、複数の磁気データq1~qNの示す座標と、球面S2との誤差を正確に表現するように定められる。従って、歪評価行列Eの各成分は、最大固有値λE1に対応する固有ベクトルuE1の向きと、球面S2からの誤差が大きな磁気データが多く存在する領域を基準点wKEの示す座標から表したベクトルの方向とが、近づくように定められる。そして、歪評価行列Eの最大固有値λE1は、球面S2からの誤差が大きな磁気データqiについての、その誤差の大きさを表現する値となる。本実施形態は、立体SDと球面との形状の相違の程度を示す歪評価値gDEを、歪評価行列Eの最大固有値λE1の値に基づいて定める。従って、歪評価値gDEによって、球面S2からの距離が大きな座標を示す磁気データqiと、球面S2との誤差の大きさ、すなわち、立体SDと球面との形状の相違の程度を評価することが可能となる。以下に、歪評価値gDEを求める方法について説明する。まず、式87に示した関数kevは、以下の式95に変形できる。また、歪誤差ベクトルkEの第i行目の要素keは、式97によって示される6次元のベクトルke2i、及び、式98によって表される歪評価行列Eの各成分を並べた6次元のベクトルeEを用いて、以下の式96に変形できる。ここで、式99で示す行列X2を導入する。行列X2は、ベクトルke2iを転置した1行6列のベクトルと、ベクトルqi-qCを転置した1行3列のベクトルとを、各行に並べて生成された、N行9列の行列である。行列X2を用いて、式92に示した歪評価関数fSDE、cは、以下の式100に表す関数gSDeに変形される。なお、ベクトルeは、以下の式101に示すように、ベクトルeEと、以下の式102に示す3次元のベクトルeXとを並べた、9次元のベクトルである。式100に示す関数gSDeを最小化する解e=e0は、以下の式103に示す連立方程式に対して、ガウス消去法や、コレスキー分解法などを適用することで求められる。なお、式103は、式100に対して、最小二乗法を適用することによって算出される正規方程式である。このようにして得た解e0に基づき、式85の歪評価行列Eを復元する。そして、式93に示した歪評価値gDEの値、即ち、歪評価行列Eのノルムを求め、歪評価値gDEの値が歪許容値δ0以下であるか否かを判定する。なお、歪評価行列Eのノルムは、歪評価行列Eの有する3つの固有値のうち絶対値が最大となる固有値λE1の絶対値に等しいため、ヤコビ法や、冪乗法により求めることができる。歪評価値gDEが歪許容値δ0以下である場合、地磁気測定装置は、処理をステップS10のオフセット採用処理に進め、球面Sの中心点cSの座標を表すベクトルを、オフセットcOFFとして採用する。一方、歪評価値gDEが歪許容値δ0よりも大きな値である場合には、球面Sの中心点cSの座標を表すベクトルを、オフセットcOFFとして採用することはできない。この場合、地磁気測定装置は、処理をステップS6の初期楕円面生成処理に進める。このように、第1の歪判定処理は、複数の磁気データq1~qNの示す座標を近傍に有する立体SDと球面との形状の相違の程度を評価する。立体SDと球面との形状の相違が無視できる程度に小さい場合には、立体SDを球面と看做すことができるため、球面Sの中心点cSの座標を表すベクトルを、オフセットcOFFとして採用することができる。この場合には、ソフトアイアン効果が生じていないと看做すことができるため、ステップS6~S8に示す楕円補正初期楕円面生成処理、最適楕円面生成処理、及び楕円面球面変換処理を行うことなく、地磁気Bgの向きを算出することができる。すなわち、本実施形態に係る地磁気測定装置は、第1の歪判定処理を行うことにより、ソフトアイアン効果の有無を判断し、ソフトアイアン効果が生じていないと判断される場合には、楕円補正を行うことなく地磁気Bgの向きを算出する。これにより、本実施形態に係る地磁気測定装置は、地磁気Bgの向きの算出に係る計算負荷を大幅に低減することが可能となる。なお、本実施形態では、歪評価値gDEが歪許容値δ0以下である場合に、オフセット採用部610aは、球面Sの中心点cSの座標を表すベクトルをオフセットcOFFとして採用するが、球面S2の中心点cS2の座標を表すベクトルをオフセットcOFFとして採用してもよい。歪評価値gDEの値が閾値δ0以下である場合には、球面Sの中心点cSの示す座標と、球面S2の中心点cS2の示す座標とは、ほぼ等しくなり、いずれもオフセットcOFFとして採用することができるからである。なお、球面S2の中心点cS2の座標は、関数gSDeを最小化する解e0のうち式101のeXに相当する3次元のベクトルを、式102に代入した場合の、変数ベクトルcとして算出される。[11.2.第2の歪判定処理]図25を参照しつつ、ステップS9において、歪判定部900が行う第2の歪判定処理について説明する。第2の歪判定処理は、ステップS5における第1の歪判定処理における判定結果が否定である場合、すなわち、図25Aに示すように、複数の磁気データq1~qNの示す座標を近傍に有する立体SDの形状が球面とは異なる歪んだ形状であると判定された場合に、図25Bに示すように、複数の変換後磁気データs1~sNの示す座標の分布形状を評価するための処理である。すなわち、第2の歪判定処理において、歪判定部900は、複数の入力座標として、複数の変換後磁気データs1~sNの示す座標を用いる。第1の歪判定処理において、複数の磁気データq1~qNの示す座標を近傍に有する立体SDの形状が球面とは異なる歪んだ形状であると評価される場合、立体SDの歪の原因が、不均一な外部磁界Bxではなく、ソフトアイアン効果である場合がある。不均一な外部磁界Bxは存在せず、且つ、ソフトアイアン効果が生じている場合、複数の磁気データq1~qNの示す座標を近傍に有する立体SDは、楕円面VEと同一の形状を有すると看做すことができる。この場合、複数の磁気データq1~qNの示す座標を、最適楕円面補正行列TOPにより変換して得られる複数の変換後磁気データs1~sNの示す座標は、球面SEOP近傍に分布するため、複数の変換後磁気データs1~sNの示す座標を近傍に有する立体SDEの形状は、球面と看做すことができる。逆に、不均一な外部磁界Bxが存在する場合、図25Aに示すように、複数の磁気データq1~qNの示す座標を近傍に有する立体SDは、球面とは異なる歪んだ形状を有し、且つ、楕円面とも異なる形状を有する。この場合、複数の磁気データq1~qNの示す座標を、最適楕円面補正行列TOPにより変換して得られる複数の変換後磁気データs1~sNの示す座標は、図25Bに示すように球面SEOPとは異なる歪んだ形状の立体SDEの近傍に分布する。このように、第2の歪判定処理は、複数の変換後磁気データs1~sNの示す座標を近傍に有する立体SDEの形状と球面例えば、球面SEOPの形状との相違の程度を評価することにより、複数の磁気データq1~qNの示す座標を近傍に有する立体SDの形状と楕円面例えば、楕円面VEの形状との相違の程度を評価する処理である。なお、第2の歪判定処理において、立体SDEの形状が、球面の形状と等しいと看做すことができる場合、ソフトアイアン効果が生じているが、不均一な外部磁界Bxの影響は存在しないため、最適中心点cEOPの示す座標をオフセットcOFFとして採用することができる。一方、図25Bに示すように、第2の歪判定処理において、立体SDEの形状が、球面の形状とは異なる歪を持った形状であると評価される場合、複数の磁気データq1~qNは不均一な外部磁界Bxの影響を受けたものであるため、複数の磁気データq1~qNに基づいてオフセットcOFFの算出をすることはできない。以下において、第2の歪判定処理の具体的な処理内容を説明する。前述の通り、第2の歪判定処理は、歪判定処理に用いる複数の入力座標として、複数の磁気データq1~qNの示す座標を用いる代わりに、複数の変換後磁気データs1~sNの示す座標を用いる点を除いて、第1の歪判定処理と同様の処理である。つまり、第2の歪判定処理は、第1の歪判定処理で用いた複数の磁気データq1~qNの示す座標に対して、以下の式104、式105に示す複数の変換後磁気データs1~sNの示す座標の値を代入した後に、11.1節で述べた第1の歪判定処理を実行する処理である。なお、歪判定処理で用いる値のうち、9節で述べた中心点算出処理において算出される値、例えば、式83に示す立体方程式に現れる行列X等も、複数の磁気データq1~qNの示す座標代わりに、複数の変換後磁気データs1~sNの示す座標を用いて算出される。例えば、式74に示す重心qCには、以下の式106に示す、複数の変換後磁気データs1~sNの示す座標の重心sCが代入されたうえで、第2の歪判定処理が行われる。これらの値の算出は中心点算出部800において行われるが、歪判定部900が行ってもよい。また、基準点wKEには、中心点cSの示す座標の代わりに、以下の式107に示す最適中心点cEOPの示す座標を代入して、第2の歪判定処理が行われる。第2の歪判定処理は、このようにして複数の変換後磁気データs1~sNの示す座標に基づいて定められる歪評価関数fSDE、cの値を最小化して求められる歪評価行列Eから歪評価値gDEを算出することにより、立体SDEと球面との形状の相違の程度を評価する。複数の変換後磁気データs1~sNの示す座標に基づき算出された歪評価値gDEが歪許容値δ0以下である場合、ソフトアイアン効果が生じているが、不均一な外部磁界Bxの影響は存在しないため、地磁気測定装置は、処理をステップS10のオフセット採用処理に進め、最適中心点cEOPの座標を表すベクトルを、オフセットcOFFとして採用する。前述の通り、第2の歪判定処理における球面S2の中心点cS2の座標を表すベクトルをオフセットcOFFとして採用してもよい。一方、複数の変換後磁気データs1~sNの示す座標に基づき算出された歪評価値gDEが歪許容値δ0よりも大きな値である場合、不均一な外部磁界Bxの影響が存在するため、地磁気測定装置は、処理をステップS1の初期化処理に戻し、球面SEOPの最適中心点cEOPの座標を表すベクトルが、オフセットcOFFとして採用されることを防止する。なお、本実施形態において、第1の歪判定処理における歪許容値δ0と、第2の歪判定処理における歪許容値δ0とは等しい値に設定されるが、異なる値に設定してもよい。[12.第2実施形態の結論]以上に示したように、第2実施形態に係る地磁気測定装置は、歪判定部900を備えることにより、複数の磁気データq1~qNの示す座標を近傍に有する立体SDと球面との形状の相違の程度を評価する。そして、立体SDの形状が球面と看做すことができる場合、簡易な計算により地磁気Bgの向きを算出することができる。具体的には、第2実施形態に係る地磁気測定装置は、複数の磁気データq1~qNの示す座標が球面近傍に分布すると歪判定部900が判定した場合、中心点算出部800が算出する中心点cSの座標を表すベクトルを、オフセットcOFFとして採用する。そして、地磁気測定装置は、中心点cSの示す座標と、磁気データqiの示す座標とに基づいて、地磁気Bgの向きを算出する。このように、第2実施形態に係る地磁気測定装置は、3次元磁気センサ60が軟磁性材料を備えない機器1aに搭載され、ソフトアイアン効果が生じていない場合には、楕円補正を行うこと無く地磁気Bgの向きを算出するため、計算負荷を軽減することができる。また、第2実施形態に係る地磁気測定装置は、楕円面補正部200、楕円面球面変換部500、及び歪判定部900を備える。そして、複数の磁気データq1~qNの示す座標を近傍に有する立体SDの形状が球面とは異なる歪んだ形状であると歪判定部900が判定した場合に、歪判定部900は、楕円面球面変換部500の算出する複数の変換後磁気データs1~sNの示す座標を近傍に有する立体SDEと球面との形状の相違の程度を評価することで、立体SDの歪の原因が、ソフトアイアン効果によるものであるのか、または、不均一な外部磁界Bxによるものであるのか、を区別する。そして、立体SDEの形状が球面とは異なる歪んだ形状である場合、つまり、立体SDの歪の原因が不均一な外部磁界Bxである場合には、地磁気測定装置は、不均一な外部磁界Bxの影響下で測定された複数の磁気データq1~qNの示す座標に基づいてオフセットcOFFが算出されることを防止する。一方、立体SDEの形状が球面であると看做せる場合、つまり、ソフトアイアン効果が生じており、且つ、不均一な外部磁界Bxの影響がない場合、地磁気測定装置は、最適中心点cEOPの示す座標、最適楕円面補正行列TOP、及び、3次元磁気センサ60の出力する磁気データqiの示す座標に基づいて、地磁気Bgの向きを算出する。このように、楕円面補正部200、楕円面球面変換部500、及び歪判定部900は、複数の磁気データq1~qNの示す座標を近傍に有する立体SDが、球面、楕円面、及び、球面とも楕円面とも異なる歪んだ形状の立体のいずれであるかを判別する、歪形状判断部4として機能し、不正確なオフセットcOFFによる不正確な地磁気Bgの算出を防止することができる。また、第2実施形態に係る地磁気測定装置は、分布判定部700を備えることにより、複数の磁気データq1~qNの示す座標が、センサ座標系ΣSにおいて3次元的な広がりを有して分布しているか否かを判定することができる。これにより、複数の磁気データq1~qNの示す座標が2次元的または1次元的に分布している場合に、中心点算出部800が中心点cSの示す座標を算出することを防止でき、不正確な中心点cSをオフセットcOFFとして採用することを防止できる。また、複数の磁気データq1~qNの示す座標が2次元的または1次元的に分布している場合に、複数の磁気データq1~qNの示す座標の分布の形状から、楕円面の形状を特定することが困難であることが多い。従って、このような場合に、楕円面補正部200が、楕円補正を行うことを防止する。<C.変形例>本発明は上述した実施形態に限定されるものではなく、例えば、以下の変形が可能である。また、以下に示す変形例のうちの2以上の変形例を、矛盾しない範囲で適宜組み合わせることもできる。1変形例1上述した実施形態において、歪誤差ベクトルkEに用いられる基準点wKEには、式89に示した中心点cS、または式107に示した最適中心点cEOPが採用されるが、本発明はこのような形態に限定されるものでは無く、基準点wKEとして、式74に示した重心qC、または式106に示した重心sCを採用してもよい。歪評価行列Eの各成分は、基準点wKEの示す座標から見て、歪評価行列Eの最大固有値λE1に対応する固有ベクトルuE1の方向と、球面S2からの誤差が大きな磁気データqiまたは、変換後磁気データsiが多く存在する領域を示す方向とが、近づくように定められる。また、歪評価行列Eの最大固有値λE1は、基準点wKEから見て、固有ベクトルuE1の方向に存在する磁気データqiまたは、変換後磁気データsiの示す座標と球面S2との誤差の大きさを表す値となる。従って、基準点wKEをどのような値に設定しても、複数の磁気データq1~qNまたは複数の変換後磁気データs1~sNの示す座標が、基準点wKEから見て広く分布している場合には、歪評価行列Eによって、立体SDまたは立体SDEと球面S2との形状の相違の程度を評価することが可能となる。2変形例2上述した実施形態及び変形例において、第1の歪判定処理及び第2の歪判定処理は、いずれも、1つの基準点wKEを用いた歪誤差ベクトルkEに基づいて算出される歪評価値gDEにより、立体SDまたは立体SDEと球面との形状の相違の程度を評価するが、本発明はこのような形態に限定されるものではなく、2つの基準点wKEを用いて算出される異なる2つの歪誤差ベクトルkEにより、2つの歪評価値gDEを算出して、立体SDまたは立体SDEと球面との形状の相違の程度を評価してもよい。例えば、第2の歪判定処理において、まず、最適中心点cEOPを基準点wKEに採用して算出された歪評価値gDEに基づいて立体SDEと球面との形状の相違の程度を評価し、その後、重心sCを基準点wKEに採用して算出された歪評価値gDEに基づいて立体SDEと球面との形状の相違の程度を評価してもよい。この場合、2度の評価の両方において、立体SDEの形状が球面の形状と看做すことができる場合に、歪判定処理の結果を肯定としてもよい。このように、2つの基準点wKEを用いて、複数の磁気データq1~qNまたは、複数の変換後磁気データs1~sNの示す座標と球面S2との誤差の大きさを評価することにより、1つの基準点wKEのみを用いる場合に比べて、立体SDまたは、立体SDEと球面との形状の相違の程度を正確に評価することができる。3変形例3上述した実施形態及び変形例において、地磁気測定装置は、第1の歪判定処理ステップS5及び第2の歪判定処理ステップS9の双方を行うが、本発明はこのような形態に限定されるものでは無く、第1の歪判定処理または第2の歪判定処理のうち、一方のみを実施してもよい。例えば、地磁気測定装置が、第1の歪判定処理のみを行う場合、複数の磁気データq1~qNの示す座標を近傍に有する立体SDの形状が、球面と看做すことができるか否かを判定することができるため、ソフトアイアン効果が生じているか否かを判定することができる。そして、ソフトアイアン効果が生じていない場合には、楕円面補正部200による楕円補正を行うこと無く、3次元磁気センサ60が出力する磁気データqiの示す座標と、中心点算出部800の算出する球面Sの中心点cSの示す座標とにより地磁気Bgの向きを算出することができるため、計算負荷の低減が可能となる。また、例えば、地磁気測定装置が、第2の歪判定処理のみを行う場合、複数の変換後磁気データs1~sNの示す座標を近傍に有する立体SDEの形状が、球面と看做すことができるかを判定することができる。そして、第2の歪判定処理における判定結果が肯定である場合、複数の磁気データq1~qNの示す座標を近傍に有する立体SDの形状は、楕円面と看做すことができるため、3次元磁気センサ60が出力する磁気データqiの示す座標、最適楕円面補正行列TOP、及び、最適中心点cEOPに基づいて、地磁気Bgの向きを算出することができる。なお、楕円面には球面が含まれるため例えば、楕円面補正行列TEの有する3つの固有値が全て1の場合、等、ソフトアイアン効果の発生の有無にかかわらず、地磁気Bgの向きを算出することができる。一方、第2の歪判定処理における判定結果が否定である場合、複数の磁気データq1~qNは不均一な外部磁界Bxの影響を受けているため、地磁気測定装置は、オフセットcOFF及び地磁気Bgの向きの算出を防止する。4変形例4上述した実施形態及び変形例において、地磁気測定装置は、楕円面最適補正値生成部400を備え、最適楕円面補正行列TOP及び最適中心点cEOPに基づいて、磁気データqiの示す座標を、変換後磁気データsiの示す座標へと変換する楕円補正を行うが、本発明はこのような形態に限定されるものでは無く、地磁気測定装置が、楕円面最適補正値生成部400を備えずに構成されてもよい。この場合、地磁気測定装置地磁気算出部600または地磁気算出部600aは、楕円面初期補正値生成部300が生成する初期楕円面補正行列T0及び初期中心点cE0の座標を、楕円面補正行列TE及びオフセットcOFFとして採用して、楕円補正を行ってもよい。この場合、地磁気測定装置は、最適楕円面生成処理を行わないため、地磁気Bgの算出に係る計算負荷を低減することができる。なお、楕円面初期補正値生成部300は、空間Ω上の第1評価軸ξ1、第2評価軸ξ2、及び、第3評価軸ξ3の3つの異なる評価軸を用いて、3つの異なる楕円面を生成した後に、これらの3つの楕円面の中心点の相互間の距離が第1閾値Δc以下であるか否かを判定し、判定結果が肯定となる場合に限り、初期楕円面補正行列T0及び初期中心点cE0を生成する。従って、初期楕円面補正行列T0及び初期中心点cE0により表される初期楕円面VE0は、複数の磁気データq1~qNの示す座標との誤差を最小化する楕円面ではないが、複数の磁気データq1~qNの示す座標の分布形状を正確に表す。すなわち、初期楕円面補正行列T0及び初期中心点cE0を用いた楕円補正により、正確な地磁気Bgの向きを算出することができる。5変形例5上述した実施形態及び変形例において、地磁気測定装置は、第2の歪判定処理に用いられる複数の入力座標として、複数の変換後磁気データs1~sNの示す座標を適用したが、本発明はこのような形態に限定されるものでは無く、複数の変換後磁気データs1~sNの示す座標を最適中心点cEOPを起点として表した複数のベクトル、すなわち、ベクトルs1-cEOP~ベクトルsN-cEOPを複数の入力座標に適用してもよい。この場合、歪判定処理に用いられるデータのサイズを小さくすることができ、処理に必要なメモリのサイズの節約をすると共に、処理速度を向上させるという利点を有する。6変形例6上述した実施形態及び変形例において、初期楕円面生成部310は、3つの楕円面第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzの各々の係数行列及び中心点の座標を算出するものであったが、本発明はこのような形態に限定されるものでは無く、初期楕円面生成部310は、第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzのうち、2つの楕円面の各々の係数行列及び中心点の座標を算出するものであってもよい。この場合、初期楕円面生成部310は、第1楕円面生成部311、第2楕円面生成部312、及び第3楕円面生成部313のうち、少なくとも2つを備えればよい。図13において説明したように、第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzのうち2つの楕円面の形状の相違の程度具体的には、2つの楕円面の2つの中心点間の距離を評価することにより、複数の磁気データq1~qNの示す座標から楕円面の形状を特定することが困難な場合であるか否かを把握することが可能である。従って、初期楕円面生成部310が、第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzのうち、少なくとも2つの楕円面各々の係数行列及び中心点の座標を算出すれば、不適切な初期楕円面補正行列T0が生成されることを防止することができる。また、初期楕円面生成部310が2つの楕円面各々の係数行列及び中心点の座標を算出する場合、3つの楕円面各々の係数行列及び中心点の座標を算出する場合に比べて、計算負荷の軽減が可能となる。なお、初期楕円面生成部310が2つの楕円面の各々の係数行列及び中心点の座標を算出する場合、初期補正値生成部330は、当該2つの楕円面のうち、少なくとも1つの楕円面の係数行列に基づいて初期楕円面補正行列T0を算出すればよい。同様に、初期補正値生成部330は、当該2つの楕円面のうち、少なくとも1つの楕円面の中心点に基づいて、初期中心点cE0の座標を算出すればよい。7変形例7上述した実施形態及び変形例において、初期楕円面中心点判定部322は、中心点cxx、中心点cyy、及び中心点czzの3つの中心点の相互間の距離の全てが、第1閾値Δc以下であるか否か第2の条件を充足するか否かを判定するものであるが、本発明はこのような判定の方法に限定されるものでは無く、初期楕円面中心点判定部322は、中心点cxx、中心点cyy、及び中心点czzのうち2つの中心点間の距離が、第1閾値Δc以下であるか否かを判定するものであってもよい。例えば、変形例6のように、初期楕円面生成部310が、第1楕円面Vxx、第2楕円面Vyy、及び第3楕円面Vzzのうち、2つの楕円面例えは、第1楕円面Vxx及び第2楕円面Vyyの各々の中心点例えば、中心点cxx及び中心点cyyの座標を算出する場合、初期楕円面中心点判定部322は、当該2つの中心点cxx及びcyyの間の距離が、第1閾値Δc以下であるか否かを判定するものであってもよい。このような判定によっても、複数の磁気データq1~qNの示す座標の分布の形状から楕円面の形状を特定することが困難な場合であるか否かを把握することが可能であり、不適切な初期楕円面補正行列T0が生成されることを防止することができる。"