(57)【特許請求の範囲】
【請求項１】
軟磁性材料を有する部品を備えた機器に組み込まれ、互いに直交する第１軸、第２軸、及び第３軸の３方向の磁気成分であって前記軟磁性材料が生じる磁界に影響された磁気成分をそれぞれ検出し、検出結果を３軸の座標系におけるベクトルデータである磁気データとして順次出力する３次元磁気センサと、
前記３次元磁気センサが出力する複数の磁気データを取得し、取得した前記複数の磁気データに基づいて、前記３軸の座標系において、互いに形状の異なる一の楕円面及び他の楕円面の各々の中心点の座標を算出する初期楕円面生成部と、
前記一の楕円面の中心点及び他の楕円面の中心点の距離が、第１閾値以下であるか否かを判定する、初期楕円面中心点判定部と、
前記初期楕円面中心点判定部の判定結果が肯定である場合、
前記一の楕円面及び前記他の楕円面のうち少なくとも１つの楕円面を球面形状に変換する初期楕円面補正行列を算出する初期補正値生成部と、
を備えることを特徴とする地磁気測定装置。
【請求項２】
前記初期楕円面生成部は、
前記３軸の座標系において、
前記第１軸成分の二乗を表す項、前記第２軸成分の二乗を表す項、及び前記第３軸成分の二乗を表す項を含む楕円面方程式のうち、
前記第１軸成分の二乗を表す項以外の項に前記複数の磁気データの各々の座標を代入して得られる値と、前記複数の磁気データの座標のうち前記第１軸成分の二乗値と、の誤差を最小化することで、前記一の楕円面の中心点の座標を算出する第１算出部と、
前記第２軸成分の二乗を表す項以外の項に前記複数の磁気データの各々の座標を代入して得られる値と、前記複数の磁気データの座標のうち前記第２軸成分の二乗値と、の誤差を最小化することで、前記他の楕円面の中心点の座標を算出する第２算出部と、
を備えることを特徴とする、
請求項１に記載の地磁気測定装置。
【請求項３】
前記３軸の座標系において、
変数ベクトルの示す座標を起点とするベクトルであって前記磁気データの座標を終点とするベクトルを第１ベクトルとし、
前記変数ベクトルの示す座標を起点とするベクトルであって前記第１ベクトルを変数行列により変換したベクトルを第２ベクトルとし、
前記第２ベクトルの終点の座標を表すデータを変換後データとし、
前記複数の磁気データに対応する複数の前記変換後データの各々の座標と、前記変数ベクトルの示す座標を中心とする球面との誤差を表し、前記変数行列の各成分と、前記変数ベクトルの各要素とを変数とする関数を、楕円面最適化関数とし、
前記楕円面最適化関数の変数の初期値として、前記初期楕円面補正行列の各成分と、前記少なくとも１つの楕円面の中心点の座標とを適用し、その後、前記楕円面最適化関数の変数の値を逐次更新することにより、前記楕円面最適化関数を最小化する解として、最適楕円面補正行列と、最適中心点の座標とを算出する、楕円面最適補正値生成部と、
前記３次元磁気センサから出力される磁気データの示す座標を終点とし前記最適中心点の座標を起点とするベクトルを、前記最適楕円面補正行列により変換することにより、地磁気の向きを算出する、地磁気算出部と、
を更に備えることを特徴とする、
請求項１または２に記載の地磁気測定装置。
【請求項４】
前記３次元磁気センサから出力される磁気データの示す座標を終点とし前記少なくとも１つの楕円面の中心点の座標を起点とするベクトルを、前記初期楕円面補正行列により変換することにより、地磁気の向きを算出する、地磁気算出部
を更に備えることを特徴とする、
請求項１または２に記載の地磁気測定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、地磁気測定装置に関する。
【背景技術】
【０００２】
近年、携帯電話等の携帯機器や、自動車等の移動体に搭載され、地磁気を検出する３次元磁気センサが開発されている。一般に、３次元磁気センサは、磁界のベクトルを互いに直交する３方向の成分に分解してスカラー量として検知するための３つの磁気センサモジュールを備え、３つの磁気センサモジュールが各々出力するスカラー量を３つの成分とする３次元のベクトルデータを出力する。
【０００３】
３次元磁気センサが搭載される携帯電話等の機器は、着磁性を有する各種金属や、電気回路等、磁界を発生させる部品が備えられることが多い。この場合、３次元磁気センサが出力するベクトルデータは、地磁気を表すベクトルの他に、機器に搭載された部品が発する磁界等を表すベクトルも含む値となる。従って、地磁気の値を正確に知るためには、３次元磁気センサが出力するベクトルデータから、機器の部品が発する内部磁界を表すベクトルを取り除く補正処理が必要となる。このように、検出対象である地磁気の正確な値を得るために、補正処理において、３次元磁気センサから出力されるデータから取り除かれる成分をオフセットと呼ぶ。
【０００４】
内部磁界は、機器の部品が発する磁界であり、機器に対して一定の方向を向き、一定の大きさを有する。内部磁界は、機器に搭載された３次元磁気センサから見た場合には、機器をどのような姿勢に変化させた場合であっても、一定の方向と一定の大きさを有するベクトルとして表される。
一方、地磁気は、磁極北に向かう水平成分と伏角方向の鉛直成分とを有する磁界であり、地面に対して一定の方向と一定の大きさとを有する一様な磁界である。従って、地面に対して機器の姿勢を変化させる場合には、機器から見た地磁気の方向も変化することになる。すなわち、機器に搭載された３次元磁気センサから見た場合、地磁気は、機器の姿勢の変化に伴い向きを変化させる一定の大きさのベクトルとして表される。
３次元磁気センサを上下左右方向に回転させて３次元的に大きく姿勢変化させつつ複数の磁気データを取得した場合、３次元磁気センサが順次出力する複数のベクトルデータの示す複数の座標は、内部磁界を表すベクトルの示す座標を中心点とし、地磁気を表すベクトルの大きさを半径とする球面近傍に分布する。
【０００５】
特許文献１には、このような地磁気及び内部磁界の性質を利用することで、３次元磁気センサの姿勢を変化させつつ取得した複数の磁気データに基づいて、内部磁界を表す一定の向き及び大きさを有するベクトルを算出し、この内部磁界を表すベクトルをオフセットとして出力データから取り除く補正処理を行うことで、正確な地磁気の向きを算出する方法が開示されている。
【０００６】
ところで、３次元磁気センサが搭載された機器の部品が、軟磁性材料を有する場合、３次元磁気センサから順次出力されるベクトルデータの示す複数の座標は、軟磁性材料が磁化される結果として発する磁界の影響により、球面近傍には分布せず、楕円面近傍に分布する。すなわち、軟磁性材料が発する磁界の影響が無ければ、球面近傍に分布するはずの複数の座標は、軟磁性材料が発する磁界の影響により、楕円面の主軸方向に伸縮するように移動する結果、球面と同一の中心点を有する楕円面近傍に分布することになる。このように、軟磁性材料が磁化される結果として発する磁界の影響により、３次元磁気センサから順次出力されるベクトルデータの示す複数の座標が楕円面近傍に分布する現象は、ソフトアイアン効果と呼ばれている。
ソフトアイアン効果が生じている場合、楕円面近傍に存在する座標に基づいて、地磁気の正確な向きを算出することはできない。地磁気の正確な向きを算出するためには、楕円面上の座標を、球面上の座標へと移動させる座標変換、すなわち、楕円面の中心点を起点として楕円面の主軸方向に伸縮するように移動させる座標変換を行う必要がある。このような、楕円面上の座標を球面上の座標へ変換する処理を、楕円補正と呼ぶ。楕円補正を行うことで算出される座標変換後の座標から、球面の中心点の示す座標を減算することにより、地磁気の向きを算出することができる。
【０００７】
非特許文献１及び非特許文献２には、ソフトアイアン効果が生じている場合に、３次元磁気センサから出力されるベクトルデータの示す楕円面上の座標を、球面上の座標に座標変換を行うための、座標変換行列を算出する方法が開示されている。
具体的には、３次元磁気センサから順次出力される複数のベクトルデータの示す座標が楕円面上に位置することを表す連立一次方程式を定め、当該連立一次方程式の解として尤もらしい値に基づいて、座標変換行列の候補となる行列を算出する。その後、座標変換行列の候補となる行列を、座標変換後の座標と球面との誤差を表す非線形関数の値を最小化する非線形最適化演算の初期値に適用し、座標変換行列の候補となる行列の各成分を逐次更新することで、座標変換行列の最適値、すなわち、座標変換後の座標と球面との誤差が最小となる座標変換行列を算出する方法が開示されている。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開２００７−２４０２７０号公報
【非特許文献】
【０００９】
【非特許文献１】J. F. Vasconcelos, G. Elkaim, C. Silvestre, P. Oliveira, and B. Cardeira, “A Geometric Approach to Strapdown Magnetometer Calibration in Sensor Frame”, in IFAC Workshop on Navigation, Guidance and Control of Underwater Vehicles, Killaloe, Ireland, April 2008
【非特許文献２】C. C. Foster and G. H. Elkaim, “Extension of a Two-Step Calibration Methodology to Include Nonorthogonal Sensor Axes”, IEEE Transactions on Aerospace and Electronic Systems, Vol. 44, No. 3, July 2008
【発明の概要】
【発明が解決しようとする課題】
【００１０】
しかし、非特許文献１及び非特許文献２に開示された方法で算出される座標変換行列の初期値、つまり、座標変換行列の候補となる行列は、非線形関数の大局的最適解となる座標変換行列とは大きく異なる値である場合がある。
非線形最適化演算に用いられる初期値が、非線形関数の大局的最適解と大きく異なる値である場合、非線形最適化演算により求められる最適解が、大局的最適解とは異なる局所最適解となる可能性が高い。従って、非特許文献１及び非特許文献２に開示された方法で算出される座標変換行列を用いて、３次元磁気センサから出力されるベクトルデータを座標変換しても、正確な地磁気の向きを求めることはできない可能性が高いという問題点があった。
【００１１】
本発明は、上述した事情を鑑み、非線形最適化問題の対象となる非線形関数の大局的最適解に近い正確な初期値を算出し、当該初期値に基づいて座標変換行列を算出することにより、正確な地磁気の方向を算出することを解決課題とする。
【課題を解決するための手段】
【００１２】
以下、本発明について説明する。なお、本発明の理解を容易にするために本実施形態、変形例、及び添付図面の参照符号を括弧書きにて付記するが、それにより本発明が本実施形態に限定されるものではない。
【００１３】
上述した課題を解決するため、本発明に係る地磁気測定装置は、軟磁性材料を有する部品を備えた機器に組み込まれ、互いに直交する第１軸、第２軸、及び第３軸の３方向の磁気成分をそれぞれ検出する３次元磁気センサと、前記３次元磁気センサから、３軸の座標系において表現される３次元のベクトルデータとして順次出力される磁気データ(ｑｉ)を蓄積する蓄積部と、複数の前記磁気データ(ｑ１〜ｑＮ)の各々が示す３軸の座標を近傍に有し、互いに形状の異なる第１楕円面(Ｖｘｘ)、第２楕円面(Ｖｙｙ)、及び第３楕円面(Ｖｚｚ)のうち、少なくとも２つの楕円面の各々の中心点を表す３軸の座標を算出する初期楕円面生成部(310)と、前記少なくとも２つの楕円面の各々の中心点を表す３軸の座標の相互間の距離が、第１閾値(Δｃ)以下であるか否かを判定する、初期楕円面中心点判定部(322)と、前記初期楕円面中心点判定部の判定結果が肯定である場合、前記少なくとも２つの楕円面のうち、少なくとも１つの楕円面の形状を表す係数行列に基づいて、初期楕円面補正行列(Ｔ０)を算出すると共に、前記少なくとも１つの楕円面の中心点を表す３軸の座標に基づいて、初期中心点(ｃＥ０)の３軸の座標を算出する、初期補正値生成部(330)と、を備えることを特徴とする。
【００１４】
ソフトアイアン効果が生じている場合、複数の磁気データの示す座標は、楕円面近傍に分布するため、地磁気の向きを算出するためには、楕円面近傍に分布する座標を、当該楕円面と同一の中心点を有する球面近傍に分布する座標へと変換する座標変換（つまり、楕円補正）を行う必要がある。このような座標変換を行うための行列を算出するためには、複数の磁気データを近傍に有し、複数の磁気データの分布形状を正確に表現する楕円面の形状を特定することが必要となる。
【００１５】
しかし、複数の磁気データの示す座標の分布の形状から、楕円面の形状を特定することが困難な場合、例えば、複数の磁気データの示す座標の分布が広がりを有さず特定の座標の近傍に分布する場合であっても、複数の磁気データの示す座標を近傍に有するような楕円面を強引に算出することができる。
このような楕円面は、複数の磁気データの分布形状を正確に表現したものではない。複数の磁気データの示す座標の分布形状を正確に表さない不適切な楕円面に基づいて生成される座標変換行列により楕円補正を行う場合、正確な地磁気の向きを算出することは困難である。従って、複数の磁気データの示す座標の分布の形状から、楕円面の形状を特定することが困難な場合には、楕円補正を行うための座標変換行列が算出されることを防止する必要がある。
【００１６】
この発明によれば、初期楕円面生成部は、複数の磁気データの示す座標を近傍に有する楕円面として、第１楕円面、第２楕円面、及び、第３楕円面の３つの異なる楕円面のうち、少なくとも２つの楕円面の各々の中心点の座標を算出する。そして、初期楕円面中心点判定部は、２つの楕円面各々の中心点の相互間の距離が、第１閾値以下であることを判定する。
初期楕円面中心点判定部の判定結果が肯定である場合、初期楕円面生成部が生成する少なくとも２つの楕円面は、いずれも、複数の磁気データの示す座標を近傍に有し、初期楕円面生成部が生成する少なくとも２つの楕円面各々の中心点の示す座標は、同一と看做すことができる。従って、初期楕円面中心点判定部の判定結果が肯定である場合、初期楕円面生成部が生成する少なくとも２つの楕円面は、同一の形状を有すると看做すことができる。
【００１７】
複数の磁気データの示す座標が、楕円面の形状を特定することができるような形状に分布している場合、初期楕円面生成部において生成される少なくとも２つの異なる楕円面は、複数の磁気データの分布により特定される楕円面の形状と、同一と看做すことのできる形状を有する楕円面として、各々算出される。
逆に、複数の磁気データの示す座標の分布形状から、楕円面の形状を特定することが困難な場合、初期楕円面生成部が生成する少なくとも２つの楕円面の形状の各々は、複数の磁気データの示す座標を近傍に有するという条件のみに基づいて定められる。この場合、少なくとも２つの楕円面の形状は、同一の形状と看做すことができない可能性が高く、少なくとも２つの楕円面の中心点の各々が示す座標も、同一の座標とは看做すことができない可能性が高い。
【００１８】
本発明に係る地磁気測定装置は、第１楕円面、第２楕円面、及び第３楕円面のうち、少なくとも２つの楕円面が同一の形状と看做すことができることを判定したうえで、初期楕円面補正行列を生成する。従って、複数の磁気データの示す座標の分布の形状から、楕円面の形状を特定することが困難な場合に、不適切な初期楕円面補正行列が生成されることを防止することができる。
【００１９】
また、上述した地磁気測定装置において、前記初期楕円面生成部(310)は、複数の前記磁気データの各々が示す３軸の座標が、楕円面の近傍に確率的に分布すると仮定して、前記楕円面を、前記第１軸成分の二乗を表す項(ｘ２)、前記第２軸成分の二乗を表す項(ｙ２)、及び前記第３軸成分の二乗を表す項(ｚ２)を含む楕円面方程式により表現したとき、複数の前記磁気データの各々が表す３軸の座標を前記楕円面方程式のうち前記第１軸成分の二乗を表す項以外の項に代入して得られる値と、複数の前記磁気データの各々が表す３軸の座標のうち前記第１軸成分の二乗値との誤差を最小化する前記第１楕円面の中心点(ｃｘｘ)を表す３軸の座標を算出する、第１楕円面生成部(311)、複数の前記磁気データの各々が表す３軸の座標を前記楕円面方程式のうち前記第２軸成分の二乗を表す項以外の項に代入して得られる値と、複数の前記磁気データの各々が表す３軸の座標のうち前記第２軸成分の二乗値との誤差を最小化する前記第２楕円面の中心点(ｃｙｙ)を表す３軸の座標を算出する、第２楕円面生成部(312)、及び、複数の前記磁気データの各々が表す３軸の座標を前記楕円面方程式のうち前記第３軸成分の二乗を表す項以外の項に代入して得られる値と、複数の前記磁気データの各々が表す３軸の座標のうち前記第３軸成分の二乗値との誤差を最小化する前記第３楕円面の中心点(ｃｚｚ)を表す３軸の座標を算出する、第３楕円面生成部(313)のうち、少なくとも２つを備えることを特徴とする、ことが好ましい。
【００２０】
複数の磁気データの示す座標が、楕円面の形状を特定することができるような形状に分布している場合、複数の磁気データの示す座標と楕円面との誤差をどのような形式により表現しても、複数の磁気データの示す座標と楕円面との誤差を最小化することにより生成される楕円面は、複数の磁気データの分布により特定される楕円面の形状と同一と看做すことのできる形状に定められる。
逆に、複数の磁気データの示す座標の分布形状から、楕円面の形状を特定することが困難な場合、複数の磁気データの示す座標と楕円面との誤差を最小化することにより生成される楕円面の形状は、誤差の表現形式に依存した形状となる。
【００２１】
この発明によれば、複数の磁気データの示す座標と楕円面との誤差を第１軸成分の二乗値に基づいて表現したときに当該誤差を最小化する第１楕円面、複数の磁気データの示す座標と楕円面との誤差を第２軸成分の二乗値に基づいて表現したときに当該誤差を最小化する第２楕円面、及び、複数の磁気データの示す座標と楕円面との誤差を第３軸成分の二乗値に基づいて表現したときに当該誤差を最小化する第３楕円面のうち、少なくとも２つの楕円面を生成する。すなわち、初期楕円面生成部の生成する少なくとも２つの楕円面の各々は、異なる形式で表現された誤差を各々最小化することで定められる。
【００２２】
従って、複数の磁気データの示す座標の分布形状から楕円面の形状を特定することが可能な場合、初期楕円面生成部の生成する少なくとも２つの楕円面は、複数の磁気データの示す座標の分布形状により特定される楕円面の形状と同一と看做すことのできる形状に定められる。
一方、複数の磁気データの示す座標の分布形状から、楕円面の形状を特定することが困難な場合、初期楕円面生成部の生成する少なくとも２つ楕円面は、誤差の表現形式に依存した異なる形状を有する。この場合、少なくとも２つ楕円面の各々の中心点は、同一の座標と看做すことはできず、初期楕円面補正行列は生成されない。
このように、本発明によれば、複数の磁気データの示す座標の分布の形状から、楕円面の形状を特定することが困難な場合に、不適切な初期楕円面補正行列が生成されることを防止することができる。
【００２３】
また、上述した地磁気測定装置は、３次元の変数ベクトル(ｃ)の示す３軸の座標を起点として前記磁気データ(ｑｉ)の示す３軸の座標を表す３次元のベクトルを第１変数ベクトル(ｑｉ−ｃ)とし、前記第１変数ベクトルを、変数行列(Ｔ)により変換して得られる３次元のベクトルを第２変数ベクトル(ｓＸ−ｃ)とし、前記第２変数ベクトルの示す３軸の座標を表すデータを変換後データ(ｓＸｉ)とし、複数の前記変換後データ(ｓＸ１〜ｓＸＮ)の各々が示す３軸の座標と、前記変数ベクトルの示す３軸の座標を中心とする球面との誤差を表し、前記変数行列の各成分と、前記変数ベクトルの各要素とを変数とする関数を、楕円面最適化関数(ｆＥＬ)としたとき、前記楕円面最適化関数の変数の初期値として、前記初期楕円面補正行列の各成分と、前記初期中心点の示す３軸の座標とを適用し、その後、前記楕円面最適化関数の変数の値を逐次更新することにより、前記楕円面最適化関数を最小化する解として、最適楕円面補正行列(ＴＯＰ)と、最適中心点(ｃＥＯＰ)の示す３軸の座標とを算出する、楕円面最適補正値生成部(400)と、前記３次元磁気センサから出力される磁気データ(ｑｉ)の示す３軸の座標を、前記最適中心点の示す３軸の座標から表した３次元のベクトル(ｑｉ−ｃＥＯＰ)を、前記最適楕円面補正行列により変換することにより(ｓｉ−ｃＥＯＰ)、地磁気(Ｂｇ)の向きを算出する、地磁気算出部(600)と、を更に備えることを特徴とすることが好ましい。
【００２４】
この発明によれば、楕円面最適化関数の変数の初期値として、初期楕円面補正行列の各成分と、初期中心点の示す座標とを適用する。
前述の通り、初期楕円面補正行列及び初期中心点は、初期楕円面生成部が生成する少なくとも２つ楕円面に基づいて生成される値であり、複数の磁気データの示す座標の分布形状を正確に表現する楕円面に基づいて定められる値である。
一方、楕円面最適化関数を最小化する非線形最適化演算は、複数の磁気データの示す座標との誤差を最小化する楕円面を算出するための演算である。つまり、楕円面最適化関数の大局的最適解は、複数の磁気データの示す座標の分布形状を最も正確に表現する楕円面の形状を表す行列及び当該楕円面の中心点となる。従って、初期楕円面補正行列及び初期中心点は、楕円面最適化関数の大局的最適化に近い値を有する正確な値である。
すなわち、本発明に係る地磁気測定装置は、大局的最適解に近い正確な値を、非線形最適化演算の初期値として適用することにより、非線形最適化演算が、局所最適解を算出することを防止し、正確な地磁気の向きを算出することができる。
【００２５】
また、上述した地磁気測定装置において、前記３次元磁気センサから出力される磁気データ(ｑｉ)の示す３軸の座標を、前記初期中心点(ｃＥ０)の示す３軸の座標から表した３次元のベクトル(ｑｉ−ｃＥ０)を、前記初期楕円面補正行列(Ｔ０)により変換することにより(ｓｉ−ｃＥ０)、地磁気(Ｂｇ)の向きを算出する、地磁気算出部(600)を更に備えることを特徴とすることが好ましい。
【００２６】
この発明によれば、簡易な計算により楕円補正を行うことができ、地磁気の向きの算出に係る計算負荷を軽減することが可能となる。
【図面の簡単な説明】
【００２７】
【図１】本発明の実施形態に係る３次元磁気センサが測定する磁界の概要を説明する概念図である。
【図２】本発明の実施形態に係る３次元磁気センサが測定する地磁気と内部磁界とについて説明する概念図である。
【図３】本発明の実施形態に係る３次元磁気センサが測定する地磁気と内部磁界とについて説明する概念図である。
【図４】本発明の実施形態に係る３次元磁気センサが測定する着磁磁界について説明する概念図である。
【図５】本発明の実施形態に係る３次元磁気センサが測定する着磁磁界について説明する概念図である。
【図６】本発明の実施形態に係る３次元磁気センサが測定する着磁磁界について説明する概念図である。
【図７】本発明の実施形態に楕円面補正行列について説明する概念図である。
【図８】本発明の実施形態に係る３次元磁気センサが搭載される機器の構成を示すブロック図である。
【図９】本発明の実施形態に係る地磁気測定装置の構成を示す機能ブロック図である。
【図１０】本発明の実施形態に係る楕円面初期補正値生成部の構成を示す機能ブロック図である。
【図１１】本発明の実施形態に係る第１楕円面、第２楕円面、及び第３楕円面について説明する概念図である。
【図１２】本発明の実施形態に係る初期楕円面補正行列について説明する概念図である。
【図１３】本発明の実施形態に係る第２の条件について説明する概念図である。
【図１４】楕円補正において回転が発生する場合について説明する概念図である。
【図１５】第２実施形態に係る地磁気測定装置の構成を示す機能ブロック図である。
【図１６】第２実施形態に係る３次元磁気センサが測定する地磁気、内部磁界、着磁磁界、及び外部磁界について説明する概念図である。
【図１７】第２実施形態に係る３次元磁気センサが測定する外部磁界について説明する概念図である。
【図１８】第２実施形態に係る３次元磁気センサが測定する外部磁界について説明する概念図である。
【図１９】第２実施形態に係る地磁気測定装置の動作を表すフローチャートある。
【図２０】第２実施形態に係る中心点算出処理を説明する概念図である。
【図２１】第２実施形態に係る磁気データ分布判定処理を説明する概念図である。
【図２２】第２実施形態に係る磁気データ分布判定処理を説明する概念図である。
【図２３】第２実施形態に係る歪判定処理を説明する概念図である。
【図２４】第２実施形態に係る歪判定処理を説明する概念図である。
【図２５】第２実施形態に係る歪判定処理を説明する概念図である。
【発明を実施するための形態】
【００２８】
＜Ａ．第１実施形態＞
以下、本発明の実施の形態を説明する。
【００２９】
［１．  ３次元磁気センサが検出する磁界の概要］                        
本実施形態では、３次元磁気センサが検出する磁界として、検出対象である地磁気の他に、３次元磁気センサが搭載される機器に備えられる部品が発する磁界である内部磁界、及び部品の備える軟磁性材料が機器外部からの磁界により磁化される結果生じる着磁磁界の存在を想定する。
以下、図１乃至図５を用いて、本実施形態において想定するこれら３種類の磁界の概要、及び、これらの磁界を３次元磁気センサが検出した場合に３次元磁気センサから出力されるベクトルデータについて説明する。
【００３０】
図１は、測定対象である地磁気Ｂｇ、３次元磁気センサの搭載された機器１が備える部品２が発する内部磁界Ｂｉ、及び、部品２のうち軟磁性材料２１が発する着磁磁界Ｂｍについて説明する図である。
【００３１】
地磁気Ｂｇは、磁極北に向かう一様な向き及び大きさを有する磁界である。地磁気Ｂｇの向き及び大きさは、厳密には地域によって異なるが、例えば、異なる都市に移動する等の大きな移動をしない場合には、一様な向き及び大きさを有する。
内部磁界Ｂｉは、機器１に備えられた部品２の発する磁界であり、機器１から見て一定の向き及び一定の大きさを有する磁界である。すなわち、機器１の姿勢がどのように変化しても、内部磁界Ｂｉは、３次元磁気センサ６０によって、一定の向き及び大きさを有する磁界として検知される。
着磁磁界Ｂｍは、機器１の外部の物体から発せられる磁界（すなわち、地磁気Ｂｇ）の影響により軟磁性材料２１が磁化される結果、軟磁性材料２１が発する磁界である。従って、着磁磁界Ｂｍは、地磁気Ｂｇの向き及び大きさや、軟磁性材料２１の材質、大きさ、形状等に依存して、向き及び大きさを変化させる磁界である。
【００３２】
説明の便宜上、図１に示すような地上座標系ΣＧ及びセンサ座標系ΣＳを導入する。図１に記載された各ベクトルの左上に付された添字Ｇは、当該ベクトルが地上座標系ΣＧにおいて表現されたベクトルであることを意味する。
地上座標系ΣＧは、地上に固定された座標系であり、地上の任意の一点を原点とし、互いに直交する３つの方向、例えば、東、北、及び鉛直上向きを、それぞれｘ軸、ｙ軸、及びｚ軸とする座標系である。
センサ座標系ΣＳは３次元磁気センサ６０に固定された座標系であり、３次元磁気センサ６０の有する３つの各々のセンサモジュールの出力する値を、それぞれ、ｘ軸（第１軸）、ｙ軸（第２軸）、ｚ軸（第３軸）上にプロットするように設けられた座標系である。すなわち、３次元磁気センサ６０が出力する磁気データは、センサ座標系ΣＳのベクトルデータとして表現される。なお、図１に示す姿勢μは、地上座標系ΣＧにおけるセンサ座標系ΣＳの各軸の向き（すなわち、地上座標系ΣＧにおける、３次元磁気センサ６０の向き）を表す。
以下において、姿勢μを変化させた場合、地上座標系ΣＧ及びセンサ座標系ΣＳの各々において、内部磁界Ｂｉ及び着磁磁界Ｂｍの向きがどのように変化するかについて説明する。
【００３３】
まず、図２及び図３を参照しながら、地上座標系ΣＧ及びセンサ座標系ΣＳにおいて、内部磁界Ｂｉ及び地磁気Ｂｇがどのように表されるかについて説明する。なお、図２及び図３においては、簡単のために、機器１は軟磁性材料２１を備えず、着磁磁界Ｂｍが存在しない場合を想定する。
図２は、内部磁界Ｂｉ及び地磁気Ｂｇの向き及び大きさを、地上座標系ΣＧにおいて表した図である。機器１の姿勢μが、姿勢μ１から姿勢μ２に変化した場合、内部磁界ＧＢｉは、大きさは一定であるが、向きは姿勢μの変化に伴い変化する。一方、地磁気ＧＢｇの向き及び大きさは一定である。
【００３４】
図３は、内部磁界Ｂｉ及び地磁気Ｂｇの向き及び大きさを、センサ座標系ΣＳにおいて表現した図である。具体的には、図３は、機器１の姿勢μをμ１〜μＮと変化させつつ磁界を測定した場合に、３次元磁気センサ６０が出力するＮ個の磁気データｑ１〜ｑＮの示す座標を、センサ座標系ΣＳにおいてプロットした図である（Ｎは精度よくオフセットを導出するために必要な磁気データの規定測定回数を表す９以上の自然数）。ここで、図３に記載された各ベクトルの左上に付された添字Ｓは、当該ベクトルがセンサ座標系ΣＳにおいて表されたベクトルであることを意味する。
【００３５】
センサ座標系ΣＳにおいて、内部磁界Ｂｉは、一定の向き及び一定の大きさを持つベクトルＳＢｉ（センサ座標系ΣＳの原点から中心点ｃＯＧへと向かうベクトル）として表現される。一方、地磁気Ｂｇは、大きさは一定であるが、向きは３次元磁気センサ６０の姿勢μに伴い変化する。つまり、地磁気Ｂｇは、機器１の姿勢μに依存した向きの、一定の大きさのベクトルＳＢｇ（μ）として表現される。従って、ベクトルＳＢｇ（μ）の起点を中心点ｃＯＧに配置して姿勢μを変化させた場合、ベクトルＳＢｇ（μ）の終点は、中心点ｃＯＧを中心点とし、地磁気Ｂｇの大きさを半径とする、球面ＳＧ上の座標を示す。
【００３６】
センサ座標系ΣＳにおいて、複数の磁気データｑ１〜ｑＮの示す座標の各々は、内部磁界ＳＢｉ及び地磁気ＳＢｇの和を表すため、複数の磁気データｑ１〜ｑＮの示す座標は、球面ＳＧ上に分布する。なお、３次元磁気センサ６０の測定値は、測定誤差を有するため、複数の磁気データｑ１〜ｑＮの示す座標は、厳密には、球面ＳＧの近傍に確率的に分布する。
従って、磁気データｑｉの示す座標から内部磁界ＳＢｉを引き算することにより、センサ座標系ΣＳにおける地磁気ＳＢｇの向き及び大きさを算出することが可能となる。
このように、検出対象である地磁気Ｂｇの正確な向きを得るために、磁気データｑｉの示す座標から、３次元磁気センサ６０の出力する内部磁界Ｂｉを表す球面ＳＧの中心点ｃＯＧの示す座標を減算する処理を、補正処理と呼ぶ。
また、補正処理において磁気データｑｉから取り除かれるベクトルをオフセットｃＯＦＦと呼ぶ。すなわち、オフセットｃＯＦＦは、内部磁界を表すベクトルＳＢｉであり、センサ座標系ΣＳにおいて、原点から球面ＳＧの中心点ｃＯＧを示すベクトルとして表される。
【００３７】
機器１が軟磁性材料２１を備える場合、地磁気Ｂｇの影響を受けて軟磁性材料２１が磁化される結果、軟磁性材料２１は着磁磁界Ｂｍを発する。機器１の姿勢μを変化させた場合、着磁磁界Ｂｍの向き及び大きさが、地上座標系ΣＧ及びセンサ座標系ΣＳの各々においてどのように変化するかについて、図４及び図５を参照しながら説明する。
【００３８】
図４は、地上座標系ΣＧにおける、着磁磁界Ｂｍの向き及び大きさ例示した図である。図４は、機器１が、センサ座標系ΣＳのｘ軸に平行な長辺２１１ａ及び２１１ｂと、センサ座標系ΣＳのｙ軸に平行な短辺２１２ａ及び２１２ｂとを有する直方体の軟磁性材料２１を備え、軟磁性材料２１が、センサ座標系ΣＳのｘ軸上に位置するように配置されている場合を例示している。
【００３９】
着磁磁界Ｂｍは、軟磁性材料２１が地磁気Ｂｇにより磁化された結果として生じる磁界であり、機器１の姿勢μ及び軟磁性材料２１の材質、大きさ、形状等に依存して、向き及び大きさを変化させる磁界である。機器１の姿勢μが姿勢μ１から姿勢μ２に変化した場合、着磁磁界ＧＢｍの向き及び大きさは、ＧＢｍ（μ１）からＧＢｍ（μ２）へと変化する。例えば、機器１の姿勢μが姿勢μ１の場合、軟磁性材料２１は、軟磁性材料２１の短辺２１２ａから他方の短辺２１２ｂへと向かう着磁磁界ＧＢｍ（μ１）を発生させ、機器１の姿勢μが姿勢μ２の場合、軟磁性材料２１は、軟磁性材料２１の長辺２１１ａから他方の長辺２１１ｂへと向かう着磁磁界ＧＢｍ（μ２）を発生させる。
３次元磁気センサ６０が検出する着磁磁界ＧＢｍ（μ）の向き及び大きさは、姿勢μと、センサ座標系ΣＳにおける軟磁性材料２１の位置ＳＰｍとに依存する。例えば、図４の場合、３次元磁気センサ６０は、着磁磁界ＧＢｍ（μ１）を、地磁気ＧＢｇ（μ１）と同じ方向を向いた磁界として測定する。また、３次元磁気センサ６０は、着磁磁界ＧＢｍ（μ２）を、地磁気ＧＢｇ（μ２）とは逆の方向を向いた磁界として測定する。
なお、軟磁性材料２１は、地磁気Ｂｇの他に、部品２が発するセンサ座標系ΣＳから見て一定の向き及び大きさの磁界により磁化される。センサ座標系ΣＳから見て一定の向き及び大きさの磁界により磁化される結果として軟磁性材料２１が発する磁界は、機器１の姿勢μが変化した場合であっても一定の方向及び大きさを有する。このような、軟磁性材料２１が磁化された結果発する磁界のうち、姿勢μが変化しても一定の向き及び大きさを有する磁界は、前述した内部磁界Ｂｉに包含されるものとする。
【００４０】
図５は、機器１が姿勢μ１となるときに測定される磁気データｑ１と、姿勢μ２となるときに測定される磁気データｑ２とを、センサ座標系ΣＳにおいてプロットした図である。
磁気データｑ１は、中心点ｃＯＧの示す座標を起点とし、地磁気ＳＢｇ（μ１）と同じ向きを有する着磁磁界ＳＢｍ（μ１）と、地磁気ＳＢｇ（μ１）とを加算したベクトルＳＢＥ（μ１）により示される座標である。よって、磁気データｑ１は球面ＳＧの外側に存在する。一方、磁気データｑ２は、中心点ｃＯＧの示す座標を起点とし、地磁気ＳＢｇ（μ２）と逆向きの着磁磁界ＳＢｍ（μ２）と、地磁気ＳＢｇ（μ２）とを加算したベクトルＳＢＥ（μ２）により示される座標である。よって、磁気データｑ２は球面ＳＧの内側に存在する。
すなわち、磁気データｑ１及びｑ２は、球面ＳＧを、ベクトルＳＢｇ（μ１）方向に引き延ばし、且つ、ベクトルＳＢｇ（μ２）方向に縮めることで得られる楕円面ＶＥ上に分布する。
このように、３次元磁気センサが、軟磁性材料を備える機器に組み込まれている場合、軟磁性材料が地磁気等の機器の外部からの磁界により磁化される結果として生じる着磁磁界の影響により、３次元磁気センサの測定する複数の磁気データの示す座標が球面近傍に分布せず、楕円面近傍に分布する。このような、軟磁性材料が磁化された結果生じる磁界の影響により、複数の磁気データの示す座標が楕円面近傍に分布する現象は、ソフトアイアン効果と呼ばれている。
【００４１】
ソフトアイアン効果が生じる場合に、３次元磁気センサ６０が出力する複数の磁気データｑ１〜ｑＮの示す座標について、図６を参照しつつ説明する。
図６は、３次元磁気センサ６０の姿勢μをμ１〜μＮ（Ｎは精度よくオフセットを導出するために必要な磁気データの規定測定回数を表す９以上の自然数）と変化させつつ磁界を測定した場合に、３次元磁気センサ６０が出力するＮ個の磁気データｑ１〜ｑＮの示す座標を、センサ座標系ΣＳにおいてプロットした図である。図６は、ソフトアイアン効果により、複数の磁気データｑ１〜ｑＮの示す座標が、中心点ｃＯＧを中心とする楕円面ＶＥ上に分布する場合を想定している。なお、図６では、３次元磁気センサ６０の測定誤差を考慮していないが、測定誤差を考慮する場合、複数の磁気データｑ１〜ｑＮの示す座標は、楕円面ＶＥ上には分布せず、楕円面ＶＥの近傍に確率的に分布する。つまり、楕円面ＶＥは、複数の磁気データｑ１〜ｑＮの示す座標との誤差を最小化するようにして定められる楕円面である。
【００４２】
楕円面ＶＥの主軸を、長い順番にＬＥ１、ＬＥ２、及びＬＥ３とし、これら３つの主軸の長さをｒＥ１、ｒＥ２、及びｒＥ３とする（但し、ｒＥ１≧ｒＥ２≧ｒＥ３＞０とする）。また、球面ＳＧの半径をｒＧとする。
このとき、中心点ｃＯＧから磁気データｑｉの表す座標を示すベクトルＳＢＥ（μｉ）は、地磁気を表すベクトルＳＢｇ（μｉ）のうち、主軸ＬＥ１に平行な成分をｒＥ１／ｒＧ倍したベクトル、主軸ＬＥ２に平行な成分をｒＥ２／ｒＧ倍したベクトル、及び主軸ＬＥ３に平行な成分をｒＥ３／ｒＧ倍したベクトルの和を表すベクトルとなる。
従って、中心点ｃＯＧから磁気データｑｉの表す座標を示すベクトルＳＢＥ（μｉ）の向きと、地磁気を表すベクトルＳＢｇ（μｉ）の向きとは異なる。また、ベクトルＳＢＥ（μｉ）及びベクトルＳＢＥ（μｊ）とのなす角（つまり、２つの磁気データｑｉ及びｑｊの示す座標を、中心点ｃＯＧから見たときのなす角）と、地磁気を表すベクトルＳＢｇ（μｉ）及びベクトルＳＢｇ（μｊ）のなす角とは、異なる。このように、ソフトアイアン効果が生じている場合、磁気データｑｉの座標から中心点ｃＯＧの座標を減算しても、地磁気ＳＢｇ（μｉ）の向きを正確に求めることはできない。
【００４３】
そこで、本実施形態では、図７に示すように、楕円面ＶＥ上の座標を、半径１の球面ＳＥ上の座標へと変換する楕円面補正行列ＴＥを算出し、楕円面補正行列ＴＥによって、磁気データｑｉの示す座標を、変換後磁気データｓｉにより表される球面ＳＥ上の座標へと変換する。中心点ｃＯＧから変換後磁気データｓｉの表す座標を示すベクトルＳＢｓ（μｉ）は、地磁気を表すベクトルＳＢｇ（μｉ）と同じ方向を向く。従って、変換後磁気データｓｉの示す座標から、中心点ｃＯＧの示す座標を減算することにより、地磁気を表すベクトルＳＢｇ（μｉ）の向きを求めることができる。
このように、地磁気Ｂｇの向きを算出するために、楕円面近傍に分布する複数の磁気データの示す座標を、当該楕円面と同一の中心点を有する半径１の球面近傍に分布する複数の座標へと変換する処理を、楕円補正と呼ぶ。
【００４４】
楕円面補正行列ＴＥが行う、楕円面ＶＥ上の磁気データｑｉの示す座標から、球面ＳＥ上の変換後磁気データｓｉの示す座標への座標変換は、以下の式（１）に表される。
ここで、楕円面補正行列ＴＥは、以下の式（２）に示される３行３列の対称行列である。また、式（３）に示される３次元の変数ベクトルｑは、磁気データｑｉの座標を表すための変数ベクトルであり、式（４）に示される３次元の変数ベクトルｓは、変換後磁気データｓｉの座標を表すための変数ベクトルであり、式（５）に示される３次元の変数ベクトルｃは、中心点ｃＯＧの座標（すなわち、オフセットｃＯＦＦ）を表すための変数ベクトルである。
なお、式（１）において、ベクトル（ｑ−ｃ）は、楕円面ＶＥの中心点ｃＯＧをセンサ座標系ΣＳの原点に平行移動させて得られる楕円面上の座標を示し、ベクトル（ｓ−ｃ）は、座標系ΣＳの原点を中心とする半径１の球面上の座標を示す。
【数１】
0005866893.tif 000002
【００４５】
前述の通り、楕円面補正行列ＴＥは、楕円面ＶＥの中心点ｃＯＧを原点とする座標系において、楕円面ＶＥ上の座標を、半径１の球面ＳＥ上の座標へと変換するための行列である。つまり、楕円面補正行列ＴＥは、互いに直交する３つの固有ベクトルの各々と、楕円面ＶＥの３つの主軸の各々とが平行になり、３つの固有ベクトルに対応する３つの固有値の各々と、楕円面ＶＥの３つの主軸の各々の長さの逆数とが等しくなるように定められる。
ここで、楕円面補正行列ＴＥの有する３つの固有ベクトルをｕＴ１、ｕＴ２、及びｕＴ３とし、対応する固有値をλＴ１、λＴ２、及びλＴ３（但し、λＴ１≧λＴ２≧λＴ３＞０）とする。このとき、固有ベクトルｕＴ１は、主軸ＬＥ１と平行となり、固有ベクトルｕＴ２は、主軸ＬＥ２と平行となり、固有ベクトルｕＴ３は、主軸ＬＥ３と平行となるように定められる。また、固有値λＴ１は、主軸ＬＥ１の長さｒＥ１の逆数に等しくなり、固有値λＴ２は、主軸ＬＥ２の長さｒＥ２の逆数に等しくなり、固有値λＴ３は、主軸ＬＥ３の長さｒＥ３の逆数に等しくなるように定められる。すなわち、楕円面補正行列ＴＥは、任意のベクトルに対して、固有ベクトルｕＴ１方向の成分を固有値λＴ１倍に伸縮し、固有ベクトルｕＴ２方向の成分を固有値λＴ２倍に伸縮し、固有ベクトルｕＴ３方向の成分を固有値λＴ３倍に伸縮する行列である。
なお、楕円面補正行列ＴＥの有する３つの固有値をλＴ１、λＴ２、及びλＴ３は、いずれも正の値であり、楕円面補正行列ＴＥは、正定値行列である。
【００４６】
ところで、図７に示すように、中心点ｃＯＧから変換後磁気データｓｉにより表される座標を示すベクトルＳＢｓ（μｉ）と、地磁気を表すベクトルＳＢｇ（μｉ）との間には、ずれ角ψが生じる場合がある。この場合、変換後磁気データｓｉから、ベクトルＳＢｇ（μｉ）の正確な向きを算出することはできない。
しかし、ずれ角ψは、軟磁性材料２１と３次元磁気センサ６０との相互位置関係（すなわち、ベクトルＳＰｍの向き及び大きさ）に依存する値である。従って、ベクトルＳＰｍよりずれ角ψを特定することが可能であり、特定されたずれ角ψと、複数の変換後磁気データｓｉとから、正確な地磁気Ｂｇの向きを算出することが可能である。また、軟磁性材料２１の配置を工夫することにより、ずれ角ψを小さく抑えることも可能である。
【００４７】
以下において、楕円面ＶＥの形状を求め、楕円面補正行列ＴＥを算出することにより楕円補正を行い、正確な地磁気Ｂｇの方向を算出する方法について述べる。
【００４８】
［２．  機器構成及びソフトウェア構成］                              
図８は、本発明の第１実施形態に係る機器１の構成を示すブロック図である。
機器１は、各種の構成要素とバスを介して接続され装置全体を制御するＣＰＵ１０、ＣＰＵ１０の作業領域として機能するＲＡＭ２０、各種プログラムやデータを記憶したＲＯＭ３０、通信を実行する通信部４０、画像を表示する表示部５０、及び磁気を検出して磁気データを出力する３次元磁気センサ６０を備える。
【００４９】
３次元磁気センサ６０は、Ｘ軸地磁気センサ６１、Ｙ軸地磁気センサ６２、及びＺ軸地磁気センサ６３を備える。各センサは、ＭＩ素子（磁気インピーダンス素子）、ＭＲ素子（磁気抵抗効果素子）などを用いて構成することができる。地磁気センサＩ／Ｆ６４は、各センサから出力信号をＡＤ変換して磁気データｑを出力する。この磁気データｑは、Ｘ軸地磁気センサ６１、Ｙ軸地磁気センサ６２、及びＺ軸地磁気センサ６３からの出力の各々を、センサ座標系ΣＳの、ｘ軸、ｙ軸及びｚ軸の３成分によって示す、センサ座標系ΣＳ上のベクトルデータである。
【００５０】
ＣＰＵ１０、ＲＡＭ２０、３次元磁気センサ６０、及び磁気データ処理プログラム７０は、３次元磁気センサ６０が検出し出力する磁気データｑに基づき、正確な地磁気の向きを示す地磁気データを算出する地磁気測定装置として機能する。
表示部５０は、ＣＰＵ１０が磁気データ処理プログラム７０を実行することにより算出した地磁気の向きを、方位情報として矢印等により表示する。なお、磁気データ処理プログラム７０は、地図アプリケーション等との連携を想定するものであってもよく、表示部５０は、地磁気の向きを示す方位情報である矢印等を地図上に表示してもよい。
【００５１】
図９は、地磁気測定装置のうち、ＣＰＵ１０が磁気データ処理プログラム７０を実行することにより実現される機能を表す、機能ブロック図である。地磁気測定装置は、複数の磁気データｑ１〜ｑＮを蓄積する蓄積部１００と、最適中心点ｃＥＯＰの座標及び最適楕円面補正行列ＴＯＰを算出する楕円面補正部２００と、磁気データｑｉ、最適中心点ｃＥＯＰ、及び最適楕円面補正行列ＴＯＰに基づいて地磁気Ｂｇの方向を算出する地磁気算出部６００と、を備える。ここで、最適中心点ｃＥＯＰとは、複数の磁気データｑ１〜ｑＮの示す座標との誤差を最小化するように定められる楕円面である最適楕円面ＶＥＯＰの中心点である。また、最適楕円面補正行列ＴＯＰとは、最適楕円面ＶＥＯＰ上の座標を、最適中心点ｃＥＯＰを中心とする球面ＳＥＯＰ上の座標へと変換するための、３行３列の対称行列である。
【００５２】
蓄積部１００は、３次元磁気センサ６０から順次出力される磁気データｑ１〜ｑＮをバッファＢＵ１に蓄積する（Ｎは、精度よくオフセットを導出するために必要な磁気データの規定測定回数を表す９以上の自然数）。バッファＢＵ１はＲＡＭ２０によって形成される。
【００５３】
楕円面補正部２００は、楕円面初期補正値生成部３００と、楕円面最適補正値生成部４００とを備える。
楕円面初期補正値生成部３００は、蓄積部１００に蓄積された複数の磁気データｑ１〜ｑＮに基づいて初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０の座標を算出する。ここで、初期中心点ｃＥ０とは、蓄積部１００に蓄積された複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する初期楕円面ＶＥ０の中心点である。また、初期楕円面補正行列Ｔ０とは、初期楕円面ＶＥ０上の座標を、初期中心点ｃＥ０を中心とする球面ＳＥ０上の座標へと変換するための３行３列の対称行列である。
【００５４】
楕円面最適補正値生成部４００は、楕円面初期補正値生成部３００が出力する初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０の座標に基づいて、複数の磁気データｑ１〜ｑＮの各々が示す座標と楕円面との誤差を最小化する最適楕円面ＶＥＯＰの中心点である最適中心点ｃＥＯＰの座標と、最適楕円面ＶＥＯＰ上の座標から最適中心点ｃＥＯＰを中心点に有する球面ＳＥＯＰ上の座標への座標変換を表す最適楕円面補正行列ＴＯＰとを算出する。
複数の磁気データｑ１〜ｑＮの示す座標と最適楕円面ＶＥＯＰとの誤差が最小化され「０」となる場合、楕円面ＶＥと最適楕円面ＶＥＯＰとは一致し、最適中心点ｃＥＯＰと中心点ｃＯＧ（すなわち、内部磁界Ｂｉの示す座標）とは一致する。
【００５５】
地磁気算出部６００は、３次元磁気センサ６０から出力される磁気データｑｉの示す座標に対して、楕円面補正行列ＴＥとオフセットｃＯＦＦとを用いた楕円補正を施し、センサ座標系ΣＳにおける地磁気ＳＢｇの向き（厳密には、ベクトルＳＢｓ（μｉ）の向き）を算出する。
具体的には、地磁気算出部６００は、オフセット採用部６１０と、地磁気ベクトル計算部６２０とを備える。オフセット採用部６１０は、楕円面補正行列ＴＥとして最適楕円面補正行列ＴＯＰを採用するとともに、オフセットｃＯＦＦとして最適中心点ｃＥＯＰの座標を示すベクトルを採用する。また、地磁気ベクトル計算部６２０は、３次元磁気センサ６０から出力される磁気データｑｉに対して、楕円面補正行列ＴＥ及びオフセットｃＯＦＦを用いた楕円補正を行い、地磁気ＳＢｇの向きを算出する。
以下において、楕円面初期補正値生成部３００、楕円面最適補正値生成部４００、及び、地磁気算出部６００の詳細について、述べる。
【００５６】
［３．  初期楕円面の生成］                              
図１０は、楕円面初期補正値生成部３００の機能構成を示す機能ブロック図である。
本実施形態では、複数の磁気データｑ１〜ｑＮに基づいて初期楕円面ＶＥ０を算出する際に、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚを生成し、これら３つの楕円面に基づいて初期楕円面ＶＥ０を生成する。
以下、本実施形態における初期楕円面ＶＥ０の生成方法について具体的に説明する。
【００５７】
楕円面初期補正値生成部３００は、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚの各々の係数行列（Ｄｘｘ、Ｄｙｙ、Ｄｚｚ）及び中心点（ｃｘｘ、ｃｙｙ、ｃｚｚ）の座標を算出する初期楕円面生成部３１０、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚの各々の係数行列及び中心点に基づいて初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０の座標を算出することが適切であるか否かを判断する初期楕円面判定部３２０、及び、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚの各々の係数行列及び中心点の座標に基づいて、初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０の座標を算出する初期補正値生成部３３０を備える。
【００５８】
初期楕円面生成部３１０は、蓄積部１００に蓄積された複数の磁気データｑ１〜ｑＮに基づいて、第１楕円面Ｖｘｘの形状を表す第１楕円面係数行列Ｄｘｘと、第１楕円面Ｖｘｘの中心点ｃｘｘの座標とを算出する第１楕円面生成部３１１、複数の磁気データｑ１〜ｑＮに基づいて、第２楕円面Ｖｙｙの形状を表す第２楕円面係数行列Ｄｙｙと、第２楕円面Ｖｙｙの中心点ｃｙｙの座標とを算出する第２楕円面生成部３１２、及び、複数の磁気データｑ１〜ｑＮに基づいて、第３楕円面Ｖｚｚの形状を表す第３楕円面係数行列Ｄｚｚと、第３楕円面Ｖｚｚの中心点ｃｚｚの座標とを算出する第３楕円面生成部３１３、を備える。
以下において、第１楕円面係数行列Ｄｘｘ、第２楕円面係数行列Ｄｙｙ、及び第３楕円面係数行列Ｄｚｚと、中心点ｃｘｘの座標、中心点ｃｙｙの座標、及び中心点ｃｚｚの座標の算出方法について述べる。
【００５９】
３次元磁気センサ６０より出力される磁気データｑの示す座標を表す変数を、式（３）で表す場合、磁気データｑを表面上に有する楕円面の方程式（楕円面方程式）は、以下の式（６）で表される。なお、式（６）は楕円面を表すので、式（６）に現れる係数θｘｘ、θｙｙ、及びθｚｚはいずれも正の値となる。
【数２】
0005866893.tif 000003
【００６０】
式（６）に示される楕円面方程式は、以下の式（７）に変形される。
【数３】
0005866893.tif 000004
【００６１】
複数の磁気データｑ１〜ｑＮの示す座標の全てが、式（６）に示された楕円面上に位置する場合、以下の式（８）が成立する。
但し、ベクトルθＸＸは、式（９）に示すように、式（７）の９つの係数を並べた９次元のベクトルである。また、行列Ｒｘｘは、式（１０）に示すように、式（１３）に示す９次元のベクトルＱｘｘに対して、式（１１）に示す複数の磁気データｑ１〜ｑＮの示す座標の各々を代入して得られるＮ個のベクトルを転置したうえで各行に並べたＮ行９列の行列である。ベクトルＷｘｘは、式（１２）に示すように、複数の磁気データｑ１〜ｑＮの示す座標のうち、ｘ成分の２乗値にマイナスを掛けた値を各成分に有する、Ｎ次元のベクトルである。
【数４】
0005866893.tif 000005
【００６２】
式（８）は、ベクトルθＸＸの各要素を変数とする連立一次方程式である。従って、式（８）をベクトルθＸＸについて解くことで、式（７）の係数が決定され、複数の磁気データｑ１〜ｑＮの示す座標を表面上に有する楕円面方程式を特定できる。
しかし、３次元磁気センサ６０の測定誤差等を考慮すると、複数の磁気データｑ１〜ｑＮの示す各々の座標の全てが、式（７）に示す楕円面上に正確に一致する位置に存在することはない。従って、式（８）は解を有さず、式（８）の解としてベクトルθＸＸを算出することはできない。そこで、本実施形態では、式（８）の解として尤もらしいベクトルθＸＸを、統計的手法を用いることにより算出する。
例えば、式（７）の右辺に現れる８つの項（ｘｙ、ｘｚ、ｙ２、ｙｚ、ｚ２、ｘ、ｙ、及びｚ）をそれぞれ独立変数と看做し、式（７）の左辺に現れるｘ２を従属変数と看做したうえで、最小二乗法を用いて式（１４）に示す正規方程式を導出し、その解として、ベクトルθＸＸを求める。この正規方程式の解として表されるベクトルθＸＸは、行列（ＲｘｘＴＲｘｘ）が正則なときに、式（１５）で表すことができる。式（１５）により算出されたベクトルθＸＸを係数として式（７）に適用することにより表される楕円面を、第１楕円面Ｖｘｘと呼ぶ。
【数５】
0005866893.tif 000006
【００６３】
ここで、図１１（Ａ）に示すように、ｘｙ、ｘｚ、ｙ２、ｙｚ、ｚ２、ｘ、ｙ、及びｚの各々を変数として表す８本の軸を有する８次元の空間πｘｘと、ｘ２の値を変数として表す第１評価軸ξ１とからなる、９次元の空間Ω上に、複数の磁気データｑ１〜ｑＮの各々をプロットする。このとき、楕円面Ｖｘｘは、楕円面Ｖｘｘと、複数の磁気データｑ１〜ｑＮが空間Ω上で示す座標との、第１評価軸ξ１方向の誤差を最小化する立体（空間Ω上の８次元平面）として求められる。すなわち、楕円面Ｖｘｘは、複数の磁気データｑ１〜ｑＮを、空間πｘｘ上にプロットして得られる複数の８次元のベクトルｑπｘｘ１〜ｑπｘｘＮの各要素を式（７）の右辺に代入して得られる値ｑｘｘ１〜ｑｘｘＮと、複数の磁気データｑ１〜ｑＮの示す座標のうちｘ軸成分の２乗値ｘ１２〜ｘＮ２との誤差を最小化する立体として定められる。
【００６４】
式（７）に示される第１楕円面Ｖｘｘの方程式は、式（１７）に示す第１楕円面係数行列Ｄｘｘを用いて、式（１６）に変形される。このとき、第１楕円面Ｖｘｘの中心点ｃｘｘの座標は、式（１８）で表される。第１楕円面生成部３１１は、以上に示したように、第１楕円面係数行列Ｄｘｘ及び第１楕円面Ｖｘｘの中心点ｃｘｘを算出し、これらを出力する。なお、式（１６）が楕円面を表すための必要条件は、第１楕円面係数行列Ｄｘｘが正定値となることである。
【数６】
0005866893.tif 000007
【００６５】
次に、式（６）に示される楕円面方程式を、以下の式（１９）に変形する。
【数７】
0005866893.tif 000008
【００６６】
式（１９）に示される方程式は、ベクトルθＹＹの各要素を変数とする連立一次方程式である式（２０）に変形される。式（８）と同様に、式（２０）の解としてベクトルθＹＹを算出することはできないため、式（２０）の解として尤もらしい値として、ベクトルθＹＹを算出する。具体的には、式（２４）に示す正規方程式において、行列（ＲｙｙＴＲｙｙ）が正則なときには、ベクトルθＹＹは式（２５）により算出される。式（２５）により特定されたベクトルθＹＹを係数として式（１９）に適用することにより表される楕円面を、第２楕円面Ｖｙｙと呼ぶ。なお、ベクトルθＹＹは、式（２１）に示す９次元のベクトルであり、行列Ｒｙｙは、式（２２）に示すＮ行９列の行列であり、ベクトルＷｙｙは、式（２３）に示すＮ次元のベクトルである。
【数８】
0005866893.tif 000009
【００６７】
ここで、図１１（Ｂ）に示すように、ｘ２、ｘｙ、ｘｚ、ｙｚ、ｚ２、ｘ、ｙ、及びｚの各々を変数として表す８本の軸を有する８次元の空間πｙｙと、ｙ２の値を変数として表す第２評価軸ξ２とからなる、９次元の空間Ω上に、複数の磁気データｑ１〜ｑＮの各々をプロットする。このとき、楕円面Ｖｙｙは、楕円面Ｖｙｙと、複数の磁気データｑ１〜ｑＮが空間Ω上で示す座標との、第２評価軸ξ２方向の誤差を最小化する立体（空間Ω上の８次元平面）として求められる。すなわち、楕円面Ｖｙｙは、複数の磁気データｑ１〜ｑＮを、空間πｙｙ上にプロットして得られる複数の８次元のベクトルｑπｙｙ１〜ｑπｙｙＮの各要素を式（１９）の右辺に代入して得られる値ｑｙｙ１〜ｑｙｙＮと、複数の磁気データｑ１〜ｑＮの示す座標のうちｙ軸成分の２乗値ｙ１２〜ｙＮ２との誤差を最小化する立体として定められる。
【００６８】
式（１９）に示される第２楕円面Ｖｙｙの方程式は、式（２７）に示す第２楕円面係数行列Ｄｙｙを用いて、式（２６）に変形される。このとき、第２楕円面Ｖｙｙの中心点ｃｙｙの座標は、式（２８）で表される。第２楕円面生成部３１２は、以上に示したように、第２楕円面係数行列Ｄｙｙ及び第２楕円面Ｖｙｙの中心点ｃｙｙを算出し、これらを出力する。なお、式（２６）が楕円面を表すための必要条件は、第２楕円面係数行列Ｄｙｙが正定値となることである。
【数９】
0005866893.tif 000010
【００６９】
次に、式（６）に示される楕円面方程式を、以下の式（２９）に変形する。
【数１０】
0005866893.tif 000011
【００７０】
式（２９）に示される方程式は、ベクトルθＺＺの各要素を変数とする連立一次方程式である式（３０）に変形される。式（８）と同様に、式（２９）の解としてベクトルθＺＺを算出することはできないため、式（２９）の解として尤もらしい値として、ベクトルθＺＺを算出する。具体的には、式（３４）に示す正規方程式において、行列（ＲｚｚＴＲｚｚ）が正則なときには、ベクトルθＺＺは式（３５）により算出される。式（３５）により特定されたベクトルθＺＺを係数として式（２９）に適用することにより表される楕円面を、第３楕円面Ｖｚｚと呼ぶ。なお、ベクトルθＺＺは、式（３１）に示す９次元のベクトルであり、行列Ｒｚｚは、式（３２）に示すＮ行９列の行列であり、ベクトルＷｚｚは、式（３３）に示すＮ次元のベクトルである。
【数１１】
0005866893.tif 000012
【００７１】
ここで、図１１（Ｃ）に示すように、ｘ２、ｘｙ、ｘｚ、ｙ２、ｙｚ、ｘ、ｙ、及びｚの各々を変数として表す８本の軸を有する８次元の空間πｚｚと、ｚ２の値を変数として表す第３評価軸ξ３とからなる、９次元の空間Ω上に、複数の磁気データｑ１〜ｑＮの各々をプロットする。このとき、楕円面Ｖｚｚは、楕円面Ｖｚｚと、複数の磁気データｑ１〜ｑＮが空間Ω上で示す座標との、第３評価軸ξ３方向の誤差を最小化する立体（空間Ω上の８次元平面）として求められる。すなわち、楕円面Ｖｚｚは、複数の磁気データｑ１〜ｑＮを、空間πｚｚ上にプロットして得られる複数の８次元のベクトルｑπｚｚ１〜ｑπｚｚＮの各要素を式（２９）の右辺に代入して得られる値ｑｚｚ１〜ｑｚｚＮと、複数の磁気データｑ１〜ｑＮの示す座標のうちｚ軸成分の２乗値ｚ１２〜ｚＮ２との誤差を最小化する立体として定められる。
【００７２】
式（２９）に示される第３楕円面Ｖｚｚの方程式は、式（３７）に示す第３楕円面係数行列Ｄｚｚを用いて、式（３６）に変形される。このとき、第３楕円面Ｖｚｚの中心点ｃｚｚの座標は、式（３８）で表される。第３楕円面生成部３１３は、以上に示したように、第３楕円面係数行列Ｄｚｚ及び第３楕円面Ｖｚｚの中心点ｃｚｚを算出し、これらを出力する。なお、式（３６）が楕円面を表すための必要条件は、第３楕円面係数行列Ｄｚｚが正定値となることである。
【数１２】
0005866893.tif 000013
【００７３】
このように、初期楕円面生成部３１０は、第１楕円面係数行列Ｄｘｘ、第２楕円面係数行列Ｄｙｙ、第３楕円面係数行列Ｄｚｚ、中心点ｃｘｘの座標、中心点ｃｙｙの座標、及び中心点ｃｚｚの座標を算出し、これらを出力する。
【００７４】
初期楕円面判定部３２０は、図１０に示すように、初期楕円面係数行列判定部３２１と、初期楕円面中心点判定部３２２とを備える。初期楕円面係数行列判定部３２１及び初期楕円面中心点判定部３２２には、第１楕円面係数行列Ｄｘｘ、第２楕円面係数行列Ｄｙｙ、第３楕円面係数行列Ｄｚｚ、中心点ｃｘｘの座標、中心点ｃｙｙの座標、及び中心点ｃｚｚの座標が供給される。
初期楕円面係数行列判定部３２１は、第１楕円面係数行列Ｄｘｘ、第２楕円面係数行列Ｄｙｙ、及び第３楕円面係数行列Ｄｚｚの全てが正定値であるという条件（第１の条件）が満たされるか否かを判定する。また、初期楕円面中心点判定部３２２は、式（３９）に示すように中心点ｃｘｘと中心点ｃｙｙとの距離が第１閾値Δｃ以下であり、且つ、式（４０）に示すように中心点ｃｙｙと中心点ｃｚｚとの距離が第１閾値Δｃ以下であり、且つ、式（４１）に示すように中心点ｃｚｚと中心点ｃｘｘとの距離が第１閾値Δｃ以下であるという条件（第２の条件）が満たされるか否かを判定する。
第１の条件に係る判定結果が肯定であり、且つ、第２の条件に係る判定結果が肯定である場合、初期楕円面判定部３２０は、第１楕円面係数行列Ｄｘｘ、第２楕円面係数行列Ｄｙｙ、第３楕円面係数行列Ｄｚｚ、中心点ｃｘｘの座標、中心点ｃｙｙの座標、及び中心点ｃｚｚの座標を出力する。
一方、第１の条件または第２の条件の判定結果が否定である場合、地磁気測定装置は、処理を中断する。
【数１３】
0005866893.tif 000014
【００７５】
なお、本実施形態では、初期楕円面判定部３２０は、初期楕円面係数行列判定部３２１及び初期楕円面中心点判定部３２２を備え、第１の条件及び第２の条件の双方を満たすか否かを判定するが、本発明はこのような判定方法に限定されるものではない。
例えば、初期楕円面判定部３２０は、初期楕円面係数行列判定部３２１を備えずに構成されてもよい。この場合、初期楕円面判定部３２０は、第１の条件についての判定を行わず、第２の条件についての判定のみを行い、判定結果が肯定である場合に、第１楕円面係数行列Ｄｘｘ、第２楕円面係数行列Ｄｙｙ、第３楕円面係数行列Ｄｚｚ、中心点ｃｘｘの座標、中心点ｃｙｙの座標、及び中心点ｃｚｚの座標を出力してもよい。
【００７６】
初期補正値生成部３３０は、初期楕円面判定部３２０が行う判定の判定結果が肯定である場合、初期楕円面判定部３２０からの出力に基づいて、初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０を算出する。
ここで、初期楕円面補正行列Ｔ０とは、図１２に示すように、初期中心点ｃＥ０を中心点とする初期楕円面ＶＥ０上の座標を、初期中心点ｃＥ０を中心点とする球面ＳＥ０上の座標へと変換する行列である。初期楕円面ＶＥ０は、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び、第３楕円面Ｖｚｚのうち、少なくとも１つの楕円面に基づいて定められる楕円面であり、複数の磁気データｑ１〜ｑＮの示す座標を表面近傍に有する楕円面である。
【００７７】
初期楕円面補正行列Ｔ０の具体的な算出方法を説明するにあたり、まず、楕円面ＶＥの形状が既知であると仮定したときの、楕円面補正行列ＴＥの算出方法について説明する（段落００４３及び図７参照）。楕円面補正行列ＴＥは、楕円面ＶＥ上の座標を、楕円面ＶＥの中心点ｃＯＧを中心とする半径１の球面ＳＥ上の座標に変換する行列であり、楕円面ＶＥの形状を表す楕円面係数行列Ｄと、楕円面ＶＥの中心点ｃＯＧとに基づき算出される。
【００７８】
式（６）に示した楕円面ＶＥを表す楕円面方程式は、以下の式（４３）に示す楕円面係数行列Ｄを用いて、以下の式（４２）に変形できる。また、楕円面ＶＥの中心点ｃＯＧの示す座標は、以下の式（４４）により表される。
【数１４】
0005866893.tif 000015
【００７９】
ここで、Ｍ行Ｍ列の正定値対称行列Ｇと、Ｍ行Ｍ列の正定値対称行列Ｈとの間に、式（４５）の関係が成立する場合、行列Ｇを行列Ｈの平方根行列と呼ぶ。以下では、行列Ｈの平方根行列Ｇを、式（４６）のように、行列の２分の１乗として表現する。
【数１５】
0005866893.tif 000016
【００８０】
このとき、行列Ｈの平方根行列Ｇは、式（４７）により求められる。但し、行列Ｕ及び行列Λは、式（４８）に示すように、行列Ｈを対角化することで算出される行列である。具体的には、行列Λは、式（４９）に示すように、行列ＨのＭ個の正の固有値λＨ１〜λＨＭを対角成分とするＭ行Ｍ列の行列であり、行列Ｕは、行列Ｈの固有値λＨ１〜λＨＭに対応する固有ベクトルを正規化したうえで各列に並べたＭ行Ｍ列の回転行列である。
【数１６】
0005866893.tif 000017
【００８１】
楕円面係数行列Ｄと、楕円面補正行列ＴＥとの関係は、以上で定義した平方根行列を用いて、以下の式（５０）で表される。なお、値ｒ（Ｄ）は、以下の式（５１）に示すように、複数の磁気データｑ１〜ｑＮの示す座標を楕円面係数行列Ｄの平方根行列により変換した複数の座標と、中心点ｃＯＧとの距離の平均値を表す値である。
【数１７】
0005866893.tif 000018
【００８２】
楕円面係数行列Ｄの固有値λＤ１、λＤ２、及びλＤ３の各々は、楕円面ＶＥの主軸の長さｒＥ１、ｒＥ２、及びｒＥ３の各々の２乗値の逆数に等しい。よって、楕円面補正行列ＴＥの固有値λＴ１、λＴ２、及びλＴ３の各々は、楕円面ＶＥの主軸の長さｒＥ１、ｒＥ２、及びｒＥ３の各々の逆数に等しくなる。従って、楕円面補正行列ＴＥにより、楕円面ＶＥ上の座標を、半径１の球面ＳＥ上の座標へと変換することができる。
【００８３】
次に、初期楕円面補正行列Ｔ０、及び、初期中心点ｃＥ０の座標の算出方法について説明する。
初期楕円面補正行列Ｔ０は、第１楕円面係数行列Ｄｘｘ、第２楕円面係数行列Ｄｙｙ、及び第３楕円面係数行列Ｄｚｚのうち、少なくとも１つに基づいて算出される。また、初期中心点ｃＥ０の座標は、中心点ｃｘｘ、中心点ｃｙｙ、及び中心点ｃｚｚのうち、少なくとも１つに基づいて算出される。
前述した第２の条件を満たす場合、中心点ｃｘｘ、中心点ｃｙｙ、及び中心点ｃｚｚのうち、任意の２点間の距離は、第１閾値Δｃよりも近い。従って、第１閾値Δｃの値を十分に小さくした場合、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚは、いずれも複数の磁気データｑ１〜ｑＮの示す座標を近傍に有し、いずれの中心点も同一の座標と看做すことができるため、これら３つの楕円面は（厳密には異なる形状を有するものの）実質的に同一の楕円面と看做すことができる。この場合、初期楕円面ＶＥ０として、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚのいずれを採用してもよい。
本実施形態では、初期楕円面ＶＥ０として、第１楕円面Ｖｘｘを採用する。このとき、初期楕円面補正行列Ｔ０及び、初期中心点ｃＥ０の座標は、以下に示す式（５２）及び式（５３）により表される。
初期補正値生成部３３０は、以上のように、初期楕円面補正行列Ｔ０と初期中心点ｃＥ０の座標とを生成し、これらを出力する。
【数１８】
0005866893.tif 000019
【００８４】
なお、本実施形態では、楕円面初期補正値生成部３００は、蓄積部１００から取得した複数の磁気データｑ１〜ｑＮを、楕円面最適補正値生成部４００に対して出力するが、楕円面最適補正値生成部４００は、蓄積部１００より直接、複数の磁気データｑ１〜ｑＮを取得してもよい。
また、本実施形態では、初期楕円面ＶＥ０として、第１楕円面Ｖｘｘを採用するが、本発明はこのような形態に限定するものではない。例えば、初期楕円面ＶＥ０として、第２楕円面Ｖｙｙを採用してもよい。このとき、初期楕円面補正行列Ｔ０は式（５４）により表され、初期中心点ｃＥ０には中心点ｃｙｙが採用される。また、初期楕円面ＶＥ０として第３楕円面Ｖｚｚを採用してもよい。この場合、初期楕円面補正行列Ｔ０は式（５５）により表され、初期中心点ｃＥ０には中心点ｃｚｚが採用される。
【数１９】
0005866893.tif 000020
【００８５】
また、初期楕円面補正行列Ｔ０は、以下に示す式（５６）により算出してもよい。この場合、初期中心点ｃＥ０の座標を、以下に示す式（５７）または式（５８）により算出してもよい。
【数２０】
0005866893.tif 000021
【００８６】
ところで、初期楕円面補正行列Ｔ０を算出する方法には、以下に示す方法（以下、「対比例」と称する）も存在する（「非特許文献２」参照）。
具体的には、まず、式（６）に示す楕円面方程式を、ｘ２の項、ｙ２の項、またはｚ２の項の何れか１つの項で割り算をすることにより、楕円面方程式と等価な、以下の式（５９）に示す連立一次方程式に変形する。次に、最小二乗法を用いて、式（５９）から、以下の式（６０）に示す正規方程式を算出する。行列（ＲＴＲ）が正則なときには、楕円面の形状を示すベクトルθは、式（６０）に示す正規方程式の解として、以下の式（６１）により算出される。式（６１）により算出されるベクトルθと、式（４３）、式（４４）、及び（５０）とを用いて、初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０を算出する。
なお、例えば、式（６）に示す楕円面方程式を、ｚ２の項により割り算することで、式（５９）に示す連立一次方程式を算出する場合、ベクトルθは、以下の式（６２）に示す９次元のベクトルであり、行列Ｒは、以下の式（６３）に示す９次元のベクトルＱに対して式（１１）に示す複数の磁気データｑ１〜ｑＮの示す座標を代入して得たベクトルを転置して各行に並べて生成される以下の式（６４）に示すＮ行９列の行列であり、ベクトルＷは、以下の式（６５）に示す９次元のベクトルである。
【数２１】
0005866893.tif 000022
【００８７】
対比例は、楕円面方程式を、ｘ２の項、ｙ２の項、またはｚ２の項の何れか１つの項で割り算をすることで生成される連立一次方程式により、楕円面の形状を決定するものである。つまり、対比例は、第１評価軸ξ１、第２評価軸ξ２、及び第３評価軸ξ３のうち、１つの評価軸のみを用いて、１つの楕円面を算出することと等しい。
この場合、複数の磁気データｑ１〜ｑＮの示す座標の分布形状が同一であっても、選択される評価軸が相違すれば、算出される楕円面の形状は異なる。しかし、対比例は、１つの楕円面のみを算出するため、仮に２つ以上の評価軸を用いて２つ以上の楕円面を算出したならば把握できるはずの複数の楕円面の形状の相違について、把握することができない。例えば、２つ以上の楕円面が算出されることを前提とする、第２の条件を用いた判定を行うことができない。
従って、図１３に示す例のように、本実施形態に係る方法により、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚのうち、少なくとも２つの楕円面を算出していたならば、第２の条件についての判定結果が否定となるような場合であっても、対比例によれば初期楕円面補正行列Ｔ０が算出されることになる。
このように、対比例は、複数の磁気データｑ１〜ｑＮの示す座標の分布の形状から楕円面の形状を特定することが困難な場合であっても、複数の磁気データｑ１〜ｑＮの示す座標の分布形状を正確に表さない不適切な初期楕円面ＶＥ０に基づいた初期楕円面補正行列Ｔ０を生成する。
【００８８】
これに対して、楕円面初期補正値生成部３００は、初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０を生成するにあたり、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚの３つの楕円面を生成する。第１楕円面Ｖｘｘは、空間Ω上の第１評価軸ξ１を基準として複数の磁気データｑ１〜ｑＮの示す座標との誤差を最小化する楕円面であり、第２楕円面Ｖｙｙは、空間Ω上の第２評価軸ξ２を基準として複数の磁気データｑ１〜ｑＮの示す座標との誤差を最小化する楕円面であり、第３楕円面Ｖｚｚは、空間Ω上の第３評価軸ξ３を基準として複数の磁気データｑ１〜ｑＮの示す座標との誤差を最小化する楕円面である。すなわち、楕円面初期補正値生成部３００は、空間Ω上の３つの異なる評価軸を用いて、３つの異なる楕円面を生成する。
さらに、楕円面初期補正値生成部３００は、第１の条件及び第２の条件を用いて、算出された３つの異なる楕円面が、互いに近い形状を有することを判定する。つまり、複数の磁気データｑ１〜ｑＮの示す座標の分布の形状が楕円面の形状とは大きく異なる形状であり、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚのうち１つ以上が楕円面とは異なる形状となる場合、第１の条件は満たされない。また、図１３に示す例のように、算出された３つの楕円面各々の中心点の相互間距離が離れている場合には、第２の条件は満たされない。
このように、楕円面初期補正値生成部３００は、初期楕円面補正行列Ｔ０を生成する際に、３つの異なる楕円面を生成したうえで、第１の条件及び第２の条件が満たされているか否かを判定するため、複数の磁気データｑ１〜ｑＮの示す座標の分布の形状から、楕円面の形状を特定することが困難な場合に、複数の磁気データｑ１〜ｑＮの示す座標の分布の形状とは異なる不正確な初期楕円面ＶＥ０に基づいた、不適切な初期楕円面補正行列Ｔ０が生成されることを防止する。
【００８９】
［４．  最適楕円面の生成］          
初期楕円面補正行列Ｔ０は、式（６６）に示すように、初期中心点ｃＥ０を起点として磁気データｑｉの示す座標を終点とするベクトル（ｑｉ−ｃＥＯ）を、初期楕円面ＶＥ０の３つの主軸方向に伸縮させることで、初期中心点ｃＥ０を起点として変換後磁気データｓ０ｉの示す座標を終点とするベクトル（ｓ０ｉ−ｃＥＯ）に変換する行列である。磁気データｑｉの示す座標が、初期楕円面ＶＥ０上に存在する場合には、変換後磁気データｓ０ｉの示す座標は、初期中心点ｃＥ０を中心とする球面ＳＥ０上に位置する。
図１２に示すように、初期楕円面ＶＥ０は、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有するように定められる楕円面であるが、複数の磁気データｑ１〜ｑＮの示す座標との誤差を最小化するように定められた楕円面ではない。従って、複数の磁気データｑ１〜ｑＮの示す座標と、初期楕円面ＶＥ０との間の誤差が大きい場合、初期楕円面ＶＥ０の中心点である初期中心点ｃＥ０の示す座標と、地磁気Ｂｇを表す球面ＳＧの中心点ｃＯＧの示す座標との間の誤差も大きくなる。この場合、磁気データｑｉの示す座標を初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０により楕円補正した変換後磁気データｓ０ｉと、初期中心点ｃＥ０とに基づいて、地磁気Ｂｇの正確な向きを算出することはできない。
【数２２】
0005866893.tif 000023
【００９０】
そこで、本実施形態では、複数の磁気データｑ１〜ｑＮの示す座標との誤差を最小化する最適楕円面ＶＥＯＰに基づいて定められる最適楕円面補正行列ＴＯＰと、最適楕円面ＶＥＯＰの中心点である最適中心点ｃＥＯＰとに基づいて、磁気データｑｉの示す座標を楕円補正する。具体的には、本実施形態に係る地磁気測定装置は、楕円面補正行列ＴＥとして最適楕円面補正行列ＴＯＰを採用するとともに、オフセットｃＯＦＦとして最適中心点ｃＥＯＰの示す座標を採用することで、式（１）を、以下の式（６７）に変形し、地磁気ＳＢｇの向きを表すベクトル（ｓｉ−ｃＥＯＰ）を算出する。以下、式（６７）の右辺に現れるベクトル（ｑｉ−ｃＥＯＰ）を第１磁気ベクトルと称し、式（６７）の左辺に現れるベクトル（ｓｉ−ｃＥＯＰ）を第２磁気ベクトルと称する。
最適楕円面ＶＥＯＰは、複数の磁気データｑ１〜ｑＮの示す座標との誤差を最小化する楕円面であるため、楕円面ＶＥと、最適楕円面ＶＥＯＰとは、同一の図形と看做すことができる。よって、最適中心点ｃＥＯＰの示す座標と中心点ｃＯＧの示す座標（オフセットｃＯＦＦ）との誤差は、初期中心点ｃＥ０の示す座標と中心点ｃＯＧの示す座標との誤差に比べて小さく、最適中心点ｃＥＯＰと中心点ｃＯＧとは、同一の座標を示すと看做すことができる。このような、最適楕円面補正行列ＴＯＰと、最適中心点ｃＥＯＰとに基づいて、磁気データｑｉの示す座標を楕円補正することにより、正確な地磁気Ｂｇの向きを得ることができる。
【数２３】
0005866893.tif 000024
【００９１】
楕円面最適補正値生成部４００は、以下の式（６８）に示す楕円面最適化関数ｆＥＬ（Ｔ，ｃ）の値が最小化されるように、変数行列Ｔの各成分及び変数ベクトルｃの各要素を逐次更新する非線形最適化演算を行い、楕円面最適化関数ｆＥＬ（Ｔ，ｃ）の値が最小化されたときの変数行列Ｔ及び変数ベクトルｃを、最適楕円面補正行列ＴＯＰ及び最適中心点ｃＥＯＰとして算出する。
ここで、楕円面最適化関数ｆＥＬ（Ｔ，ｃ）は、以下の式（６９）に示す３行３列の対称行列である変数行列Ｔの各成分と、式（５）に示した変数ベクトルｃの各要素とを変数とする関数であり、以下の式（７０）として表すことができる。変数行列Ｔ及び変数ベクトルｃの初期値には、初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０が適用される。
【数２４】
0005866893.tif 000025
【００９２】
楕円面最適化関数ｆＥＬ（Ｔ，ｃ）は、式（６８）に示すように、変数ベクトルｃの示す座標を起点として複数の磁気データｑ１〜ｑＮの示す座標を表す複数の第１変数ベクトル（ｑｉ−ｃ）の各々を、変数行列Ｔにより変換することで得られる複数の第２変数ベクトルＴ（ｑｉ−ｃ）の長さの平均値が、「１」からどの程度の差異を有するかを表す関数である。
つまり、楕円面最適化関数ｆＥＬ（Ｔ，ｃ）は、複数の第２変数ベクトルＴ（ｑｉ−ｃ）を、変数ベクトルｃの示す座標が起点となるように配置した場合、複数の第２変数ベクトルの各々が示す座標と、変数ベクトルｃの示す座標を中心とする半径１の球面との誤差を表す。このとき、複数の第２変数ベクトルＴ（ｑｉ−ｃ）の示す複数の座標を表すデータを、複数の変換後データｓＸ１〜ｓＸＮと呼ぶ。楕円面最適化関数ｆＥＬ（Ｔ，ｃ）の値を最小化することにより、複数の変換後データｓＸ１〜ｓＸＮの示す座標と、変数ベクトルｃの示す座標を中心とする半径１の球面との誤差を最小化することができ、このときの複数の変換後データｓＸ１〜ｓＸＮは、複数の変換後磁気データｓ１〜ｓＮを表す。
なお、本実施形態では、説明の便宜上、第２変数ベクトルＴ（ｑｉ−ｃ）を、変数ベクトルｃの示す座標を起点として配置しているが、第２変数ベクトルＴ（ｑｉ−ｃ）を、センサ座標系ΣＳの原点が起点となるように配置しても構わない。すなわち、式（６８）は、第２変数ベクトルＴ（ｑｉ−ｃ）をセンサ座標系ΣＳの原点が起点となるように配置した場合、第２変数ベクトルＴ（ｑｉ−ｃ）の示す座標と、センサ座標系ΣＳの原点を中心とする半径１の球面との誤差を表す。また、この場合、複数の変換後磁気データｓ１〜ｓＮは、原点を中心とする半径１の球面近傍に分布する。
【００９３】
楕円面最適化関数ｆＥＬ（Ｔ，ｃ）の示す値を最小化することにより最適楕円面補正行列ＴＯＰ及び最適中心点ｃＥＯＰを算出する非線形最適化演算は、公知の方法を適宜適用すればよい。例えば、非線形最適化演算として、ニュートン法を適用してもよい。
ニュートン法、最急降下法等の、非線形最適化演算は、非線形関数の示す値を最適化（最小化または最大化）するように、非線形関数の変数の値を逐次更新する。そして、非線形最適化演算は、非線形関数の値及び変数の値が予め規定された停止規準を満たしたときに、変数の値の更新を停止し、このときの変数の値を最適解として採用する。
なお、非線形最適化演算の停止規準としては、公知の基準を適宜適用すればよい。例えば、アルミホの規準（Armijo’s  rule）を適用してもよい。
【００９４】
非線形最適化演算は、非線形関数の最適解、すなわち非線形関数を最小化（または最大化）する大局的最適解を算出するための演算である。しかし、非線形最適化演算に適用される初期値が大局的最適解とは大きく異なる値である場合、非線形最適化演算の算出する最適解が、大局的最適解とは異なる局所的最適解となることがある。大局的最適解と大きく異なる初期値を適用した場合、大局的最適解よりも初期値に近い局所的最適解が存在する可能性があり、非線形最適化演算による解の更新を繰り返すうちに、解が大局的最適解に更新される前に、局所的最適解に更新される可能性が高いからである。従って、非線形最適化演算により局所的最適解が算出されることを防止するためには、初期値として、大局的最適解にできるだけ近い値を採用する必要がある。
本実施形態は、初期楕円面ＶＥ０に基づいて、初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０を算出し、これらの値を非線形最適化演算の初期値として適用する。複数の磁気データｑ１〜ｑＮの示す座標を近傍に有するように定められた初期楕円面ＶＥ０は、複数の磁気データｑ１〜ｑＮの示す座標との誤差を最小化する楕円面ＶＥと近い形状を有する。従って、初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０の示す座標は、大局的最適解として算出されるべき値である楕円面補正行列ＴＥ及び中心点ｃＯＧの示す座標に近い値である（図７及び図１２参照）。このような初期値を用いた非線形最適化演算を行うことにより算出される最適楕円面補正行列ＴＯＰ及び最適中心点ｃＥＯＰの示す座標は、局所的最適解ではなく、大局的最適解（厳密には、大局的最適解の近似値）となる。このように、本実施形態の非線形最適化演算は、大局的最適解に近い値の適切な初期値が設定されるため、局所的最適解に陥ることなく、大局的最適解を最適解として算出することができる。
【００９５】
ところで、本実施形態では、式（６９）に示すように、非線形最適化演算において最適化の対象となる楕円面最適化関数ｆＥＬ（Ｔ，ｃ）の変数として実対称行列である変数行列Ｔを用いるが、変数行列Ｔの代わりに、対称行列に限定されない実行列一般を表す変数行列ＴＲを変数として用いて非線形最適化演算を行う方法も存在する（「非特許文献１」参照）。
しかし、対称行列に限定されない変数行列ＴＲは、任意のベクトルを変数行列ＴＲの３つの固有ベクトル方向に各々伸縮させる座標変換に加え、任意のベクトルの向きを回転させる座標変換を表す場合がある。そのため、第１変数ベクトル（ｑｉ−ｃ）を変数行列ＴＲにより変換することで得られるベクトルＴＲ（ｑｉ−ｃ）は、第２磁気ベクトルＴ（ｑｉ−ｃ）を任意角度回転させたベクトルとして算出される場合がある。すなわち、図１４に示すように、変数行列ＴＲを用いて算出される変換後磁気データＲＯｓｉの示す座標は、変数行列Ｔを用いて算出される変換後磁気データｓｉの示す座標を、球面ＳＥＯＰ上で任意角度回転させた座標として算出される。
この場合、変数行列ＴＲを用いて算出される変換後磁気データＲＯｓｉの示す座標に基づいて、地磁気Ｂｇの向きを算出することは困難である。
【００９６】
非特許文献１によれば、３次元磁気センサ６０から見て方向が既知の参照磁界を用いることにより、変数行列ＴＲが行う座標変換において生じた回転の角度を特定し、座標変換において回転が生じなかった場合の変換後磁気データｓｉの示す座標を算出している。この方法は、機器１が、参照磁界を測定する機会を有することが必要となる。
【００９７】
これに対して、本実施形態では、変数行列Ｔを実対称行列に限定している。実対称行列は、互いに直交する３つの固有ベクトルと、３つの固有ベクトルに対応する３つの固有値とを有する。そして、実対称行列によりベクトルを変換する場合、変換後のベクトルは、変換前のベクトルを実対称行列の３つの固有ベクトル方向を向く３つのベクトルの和で表したとき、３つのベクトルの各々を、その方向を変えずに対応する固有値倍だけ伸縮させた３つのベクトルの和として算出される。すなわち、実対称行列は、任意のベクトルを、当該実対称行列が有する各固有ベクトル方向に各々伸縮させる座標変換を行う行列である。
従って、実対称行列である変数行列Ｔを用いる非線形最適化演算は、最適楕円面補正行列ＴＯＰを、任意のベクトルを最適楕円面補正行列ＴＯＰの各固有ベクトル方向に各々伸縮させる座標変換を表す行列として算出するため、最適楕円面補正行列ＴＯＰは、回転を伴う座標変換は行わない。このような最適楕円面補正行列ＴＯＰにより、磁気データｑｉの示す座標を、変換後磁気データｓｉの示す座標へと変換することにより、地磁気Ｂｇの正確な向きを得ることができる。
また、変数行列ＴＲは３行３列の行列で９個の独立した成分を有し、楕円面最適化関数ｆＥＬ（ＴＲ，ｃ）の変数は１２個となるのに対して、本実施形態に係る変数行列Ｔは対称行列であるため６個の独立した成分を有し、楕円面最適化関数ｆＥＬ（Ｔ，ｃ）の変数は９個である。従って、本実施形態に係る非線形最適化演算は、変数行列ＴＲを用いる非線形最適化演算に比べて変数の数が少なく、計算負荷が低減される。
【００９８】
なお、楕円面最適補正値生成部４００は、最適楕円面補正行列ＴＯＰが正定値行列であることを判定したうえで、最適楕円面補正行列ＴＯＰ及び最適中心点ｃＥＯＰを出力してもよい。最適楕円面補正行列ＴＯＰは、楕円面上の座標を、当該楕円面の主軸方向に伸縮させることにより、球面上の座標へと変換する行列であるため、最適楕円面補正行列ＴＯＰの有する３つの固有値は、いずれも正の値となるからである。
楕円面最適補正値生成部４００が出力した最適楕円面補正行列ＴＯＰ及び最適中心点ｃＥＯＰは、蓄積部１００に格納される。
【００９９】
［５．  地磁気の算出］                              
前述の通り、地磁気算出部６００は、オフセット採用部６１０と、地磁気ベクトル計算部６２０とを備え、３次元磁気センサ６０から出力される磁気データｑｉの示す座標に対して、楕円補正を施すことにより、地磁気ＳＢｇの向きを算出する（図９参照）。以下において、地磁気算出部６００の動作を説明する。
まず、オフセット採用部６１０は、楕円面補正行列ＴＥとして最適楕円面補正行列ＴＯＰを採用するとともに、オフセットｃＯＦＦとして最適中心点ｃＥＯＰの座標を示すベクトルを採用する。これにより、本実施形態に係る地磁気測定装置は、楕円補正を表す式（１）を式（６７）に変形し、式（６７）に基づいた楕円補正を行うことが可能となる。
次に、地磁気ベクトル計算部６２０は、式（６７）に基づいて楕円補正を行い、地磁気Ｂｇの向きを算出する。具体的には、地磁気ベクトル計算部６２０は、オフセットｃＯＦＦである最適中心点ｃＥＯＰの座標を起点とし、磁気データｑｉの示す座標を終点とする第１磁気ベクトル（ｑｉ−ｃＥＯＰ）を、楕円面補正行列ＴＥである最適楕円面補正行列ＴＯＰにより変換し、第２磁気ベクトル（ｓｉ−ｃＥＯＰ）を算出する。このとき、第２磁気ベクトル（ｓｉ−ｃＥＯＰ）は、ずれ角ψを考慮しなければ、地磁気Ｂｇと同じ向きを向く。地磁気ベクトル計算部６２０は、必要な場合にはずれ角ψを考慮して（図７、及び段落００４６参照）、第２磁気ベクトル（ｓｉ−ｃＥＯＰ）から、地磁気Ｂｇの向きを算出する。
なお、本実施形態では、説明の便宜上、第２磁気ベクトル（ｓｉ−ｃＥＯＰ）を、最適中心点ｃＥＯＰの示す座標を起点として配置しているが（図７参照）、第２磁気ベクトル（ｓｉ−ｃＥＯＰ）を、センサ座標系ΣＳの原点を起点として配置しても構わない。この場合、球面ＳＥＯＰは、センサ座標系ΣＳの原点を中心とする半径１の球面を表し、変換後磁気データｓｉは原点を中心とする半径１の球面近傍に分布する。
【０１００】
［６．  第１実施形態の結論］                
以上に説明した第１実施形態では、複数の磁気データｑ１〜ｑＮの示す座標が、楕円面近傍に分布すると仮定して、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する楕円面を特定し、複数の磁気データｑ１〜ｑＮの示す座標を、当該楕円面と同一の中心点を有する球面近傍に座標変換する楕円補正を行った。
これにより、３次元磁気センサ６０を搭載する機器が軟磁性材料を備え、ソフトアイアン効果が生じている場合に、複数の磁気データｑ１〜ｑＮに基づいて、地磁気Ｂｇの正確な方向を算出することができる。
【０１０１】
なお、３次元磁気センサ６０を搭載する機器が軟磁性材料を備えず、ソフトアイアン効果が生じていない場合であっても、複数の磁気データｑ１〜ｑＮの示す座標が楕円面の近傍に分布する場合がある。例えば、３次元磁気センサ６０が備える、Ｘ軸地磁気センサ６１、Ｙ軸地磁気センサ６２、及びＺ軸地磁気センサ６３の３つセンサの感度にバラつきがある場合、センサ座標系ΣＳにおいて本来球面近傍に分布するはずの複数の磁気データｑ１〜ｑＮの示す座標は、３つのセンサの感度に応じて当該球面をセンサ座標系ΣＳの各軸方向に伸縮させた楕円面近傍に分布する。すなわち、センサ感度のバラつきがある場合、複数の磁気データｑ１〜ｑＮの示す座標は、センサ座標系ΣＳの３つの軸の向きと同じ方向を向く３つの主軸を有する楕円面近傍に分布する。
第１実施形態では、複数の磁気データｑ１〜ｑＮの示す座標が楕円面近傍に分布しているものと看做して楕円補正を行う。この楕円補正に用いる楕円面の３つの主軸は、直交関係を保ちつつ任意の方向を向くことができるため、第１実施形態に係る楕円補正は、楕円面の主軸とセンサ座標系ΣＳの３軸の方向が一致する場合についても、適用することが可能である。
従って、第１実施形態の楕円補正は、センサ感度のバラつきが原因で複数の磁気データｑ１〜ｑＮの示す座標が楕円面近傍に分布する場合、つまり、楕円面の主軸とセンサ座標系ΣＳの３軸の方向が一致する場合であっても、地磁気Ｂｇの正確な向きを算出することができる。
【０１０２】
また、第１実施形態では、初期楕円面ＶＥ０の算出に当たり、空間Ω上の第１評価軸ξ１、第２評価軸ξ２、及び第３評価軸ξ３の３つの異なる評価軸を用いて、３つの異なる楕円面を生成したうえで、センサ座標系ΣＳにおいてこれらの３つの異なる楕円面の各々の中心点の相互間の距離が全て第１閾値Δｃ以下であるか否かを判定した。そして判定結果が肯定となる場合に限り、初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０の座標を算出した。
これにより、複数の磁気データｑ１〜ｑＮの示す座標の分布の形状から、楕円面の形状を特定することが困難な場合に、複数の磁気データｑ１〜ｑＮの示す座標の分布の形状とは異なる不正確な初期楕円面ＶＥ０に基づいた、不適切な初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０の座標が生成されることを防止することが可能となり、不正確な値がオフセットとして採用されることを防止することが可能となった。
【０１０３】
また、第１実施形態では、楕円面ＶＥに近い形状を有する初期楕円面ＶＥ０に基づいて定められる初期楕円面補正行列Ｔ０の各成分及び初期中心点ｃＥ０の座標を初期値とする非線形最適化演算により、最適楕円面補正行列ＴＯＰ及び最適中心点ｃＥＯＰの座標を算出した。
非線形最適化は、初期値に、大局的最適解に近い値を採用することで、局所的最適解が最適解として算出される可能性を低減させ、大局的最適解が最適解として算出される可能性を高める。従って、本実施形態に係る非線形最適化演算は、局所的最適解が最適解として算出される可能性を低減し、不適切な最適楕円面補正行列ＴＯＰを用いた楕円補正による不正確な地磁気Ｂｇの向きが算出される可能性を低減した。
【０１０４】
また、第１実施形態では、実対称行列である変数行列Ｔを変数とする楕円面最適化関数ｆＥＬ（Ｔ，ｃ）を最小化することで、最適楕円面補正行列ＴＯＰ及び最適中心点ｃＥＯＰを算出した。これにより、最適楕円面補正行列ＴＯＰは、任意のベクトルを、最適楕円面補正行列ＴＯＰの３つの固有ベクトル方向に伸縮させる座標変換を行う行列として算出されるため、回転を伴う座標変換は行わない。
すなわち、最適楕円面ＶＥＯＰ近傍に位置する磁気データｑｉの示す座標を楕円補正することで算出される変換後磁気データｓｉの示す座標は、最適中心点ｃＥＯＰから見て、地磁気Ｂｇと同じ方向を向く座標として求められるため、変換後磁気データｓｉに基づいて、正確な地磁気Ｂｇの向きを算出することができる。
【０１０５】
＜Ｂ．第２実施形態＞                  
以下、第２実施形態について説明する。
【０１０６】
［７．  第２実施形態に係る地磁気測定装置の概要］    
第１実施形態では、３次元磁気センサ６０が測定する対象の磁界を、地磁気Ｂｇ、内部磁界Ｂｉ、及び着磁磁界Ｂｍに限定し、複数の磁気データｑ１〜ｑＮの示す座標が楕円面近傍に分布する場合を想定した。
しかし、機器１の外部に磁界を生じる物体が存在する場合、物体の発する外部磁界Ｂｘの影響により、複数の磁気データｑ１〜ｑＮの示す座標は、楕円面とは異なる形状に分布することがある。この場合には、複数の磁気データｑ１〜ｑＮの示す座標が楕円面近傍に分布すると仮定して算出された楕円面の中心点の示す座標は、３次元磁気センサ６０のオフセットを示す座標とは一致しない可能性が高く、楕円面の中心点の示す座標をオフセットとして補正を行っても正確な地磁気Ｂｇの向きを求めることはできない。
また、機器１が軟磁性材料２１を備えない場合、ソフトアイアン効果は生じないため、複数の磁気データｑ１〜ｑＮの示す座標は、図３に示したように球面近傍に分布し、楕円面近傍には分布しない。この場合、楕円補正を行わなくても、正確な地磁気Ｂｇの向きを求めることができる。
第２実施形態は、機器１の外部の物体の発する外部磁界Ｂｘが存在する場合、及び、機器１が軟磁性材料２１を備えず着磁磁界Ｂｍが存在しない場合に対応した地磁気測定装置を実現することを目的とする。
【０１０７】
図１５は、第２実施形態に係る地磁気測定装置のうち、ＣＰＵ１０が磁気データ処理プログラムを実行することにより実現される機能を表す、機能ブロック図である。第２実施形態に係る地磁気測定装置は、楕円面球面変換部５００、分布判定部７００、中心点算出部８００、及び歪判定部９００を備える点と、地磁気算出部６００の代わりに地磁気算出部６００ａを備える点とを除き、第１実施形態に係る地磁気測定装置（図９参照）と同様に構成される。
【０１０８】
楕円面球面変換部５００は、式（６７）に基づいて、最適楕円面補正行列ＴＯＰ、最適中心点ｃＥＯＰ、及び、複数の磁気データｑ１〜ｑＮから、複数の変換後磁気データｓ１〜ｓＮを算出する。具体的には、楕円面球面変換部５００は、まず、式（６７）に示ように、最適中心点ｃＥＯＰの座標を起点として磁気データｑｉの示す座標を終点とする第１磁気ベクトル（ｑｉ−ｃＥＯＰ）を、最適楕円面補正行列ＴＯＰにより変換し、最適中心点ｃＥＯＰを起点として変換後磁気データｓｉの示す座標を終点とする第２磁気ベクトル（ｓｉ−ｃＥＯＰ）を算出することで、変換後磁気データｓｉの示す座標を算出する。その後、楕円面球面変換部５００は、算出した複数の変換後磁気データｓ１〜ｓＮを蓄積部１００のバッファＢＵ２に格納する。
分布判定部７００は、センサ座標系ΣＳにおいて、複数の磁気データｑ１〜ｑＮの示す座標の分布が、３次元的な広がりを有しているか否かを判定し、判定結果を出力する。
中心点算出部８００は、センサ座標系ΣＳにおいて、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する球面Ｓの中心点ｃＳの示す座標を算出する。図３において説明したように、３次元磁気センサ６０が測定する磁界が、地磁気Ｂｇ及び内部磁界Ｂｉである場合、複数の磁気データｑ１〜ｑＮの示す座標が、球面ＳＧの近傍に分布するため、球面Ｓと球面ＳＧとは一致すると看做すことができ、球面Ｓの中心点ｃＳの示す座標が、オフセットｃＯＦＦを表す。
第２実施形態に係る地磁気測定装置は、中心点算出部８００を備えることで、着磁磁界Ｂｍが存在せず、複数の磁気データｑ１〜ｑＮの示す座標が球面近傍に分布する場合に、簡易な方法により、３次元磁気センサ６０のオフセットｃＯＦＦを算出することが可能となる。
【０１０９】
歪判定部９００は、センサ座標系ΣＳにおいて、複数の入力座標（複数の磁気データｑ１〜ｑＮの示す座標、または複数の変換後磁気データｓ１〜ｓＮの示す座標）がある立体近傍に分布すると仮定して、当該立体の形状と球面の形状との相違の程度を評価することで、立体の形状が球面と看做すことができる否かを判定し、判定結果を出力する。
外部磁界Ｂｘが存在し、複数の磁気データｑ１〜ｑＮの示す座標が球面または楕円面とは異なる歪んだ形状の立体近傍に分布する場合、３次元磁気センサ６０のオフセットｃＯＦＦの正確な値を算出することは、困難である。
第２実施形態に係る地磁気測定装置は、歪判定部９００を備えることで、外部磁界Ｂｘの影響が大きく、オフセットｃＯＦＦの算出が困難な場合に、不正確なオフセットｃＯＦＦが算出されることを防止し、不正確なオフセットを用いた補正処理により不正確な地磁気Ｂｇが算出されることを防止する。
【０１１０】
地磁気算出部６００ａは、オフセット採用部６１０の代わりにオフセット採用部６１０ａを備える点を除き、地磁気算出部６００と同様に構成される。オフセット採用部６１０ａは、オフセットｃＯＦＦとして、最適中心点ｃＥＯＰの座標を示すベクトルまたは球面Ｓの中心点ｃＳの座標を示すベクトルのいずれか一方を採用する。
そして、オフセット採用部６１０ａは、オフセットｃＯＦＦとして中心点ｃＳの座標を示すベクトルを採用する場合、楕円面補正行列ＴＥとして３行３列の単位行列Ｉを採用する。このとき、地磁気ベクトル計算部６２０は、オフセットｃＯＦＦである中心点ｃＳの示す座標と、楕円面補正行列ＴＥである単位行列Ｉとを用いて、式（１）に基づいた楕円補正を行い、地磁気Ｂｇの向きを算出する。具体的には、楕円面補正行列ＴＥは単位行列Ｉであるため、地磁気ベクトル計算部６２０は、ベクトル（ｑｉ−ｃＳ）を、地磁気ＳＢｇの向きを表すベクトルとして算出する。なお、単位行列Ｉを用いる楕円補正とは、式（１）からも明らかなとおり、単に、磁気データｑｉの示す座標からオフセットｃＯＦＦを減算する演算に過ぎず、実質的には楕円補正を行わないことを意味する。従って、楕円面補正行列ＴＥとして単位行列Ｉが採用された場合、地磁気ベクトル計算部６２０は、式（１）に基づいた演算を行わず、単純に、磁気データｑｉの示す座標からオフセットｃＯＦＦとして採用された中心点ｃＳの座標を減算する処理を行ってもよい。
一方、オフセット採用部６１０ａは、オフセットｃＯＦＦとして最適中心点ｃＥＯＰの座標を示すベクトルを採用する場合、楕円面補正行列ＴＥとして最適楕円面補正行列ＴＯＰを採用する。このとき、地磁気ベクトル計算部６２０は、オフセットｃＯＦＦである最適中心点ｃＥＯＰの示す座標と、楕円面補正行列ＴＥである最適楕円面補正行列ＴＯＰとを用いて、式（１）に基づいた楕円補正を行い、地磁気Ｂｇの向きを算出する。具体的には、地磁気ベクトル計算部６２０は、式（１）を変形した式（６７）により、第２磁気ベクトル（ｓｉ−ｃＥＯＰ）を、地磁気ＳＢｇの向きを表すベクトルとして算出する。
【０１１１】
以下では、第２実施形態に係る地磁気測定の処理の詳細を説明する前提として、外部磁界Ｂｘの性質について説明する。
【０１１２】
図１６は、３次元磁気センサ６０の検出対象である、地磁気Ｂｇ、内部磁界Ｂｉ、着磁磁界Ｂｍ、及び外部磁界Ｂｘを、地上座標系ΣＧにおいて表した概念図である。ここで、図１６に示す位置ＰＳは、地上座標系ΣＧにおけるセンサ座標系ΣＳの原点の位置（すなわち、地上座標系ΣＧにおける、３次元磁気センサ６０の位置）を表す。
第２実施形態に係る地磁気測定装置は、機器１の他に、図１６に示す機器１ａに適用することができる。ここで、機器１ａは、部品２の代わりに、軟磁性材料２１を含まない部品２ａを備える点を除き、機器１と同様の構成である。すなわち、機器１ａは、機器１と異なり、着磁磁界Ｂｍを発しない。
【０１１３】
図１６に示すように、外部磁界Ｂｘは、機器１または機器１ａの外部に存在する物体３から発せられる磁界であり、物体３との相対的位置関係によって向き及び大きさが変化する不均一な磁界である。地上座標系ΣＧにおける３次元磁気センサ６０の位置Ｐｓを変化させた場合、３次元磁気センサ６０が測定する外部磁界Ｂｘの向き及び大きさが変化する。従って、外部磁界Ｂｘは、地上座標系ΣＧにおいて、位置Ｐｓに依存して方向及び大きさの双方を変化させるベクトルＧＢｘ（Ｐｓ）として表現される。また、地上座標系ΣＧにおける３次元磁気センサ６０の姿勢μを変化させた場合、３次元磁気センサ６０が測定する外部磁界Ｂｘの向きが変化する。
【０１１４】
図１７は、３次元磁気センサ６０の位置ＰｓをＰＳ１〜ＰＳＮと変化させると共に、姿勢μをμ１〜μＮと変化させて磁界を測定したときの、３次元磁気センサ６０が出力する複数の磁気データｑ１〜ｑＮを、センサ座標系ΣＳにおいてプロットした図である。
なお、図１７では、簡単のために、着磁磁界Ｂｍは存在せず、内部磁界Ｂｉ、地磁気Ｂｇ、及び外部磁界Ｂｘのみが存在する場合を想定する。
【０１１５】
外部磁界Ｂｘは、センサ座標系ΣＳにおいて、３次元磁気センサ６０の位置Ｐｓに依存して方向及び大きさの双方を変化させ、３次元磁気センサ６０の姿勢μに依存して方向を変化させるベクトルＳＢｘ（μ、Ｐｓ）として表現される。
３次元磁気センサ６０が、内部磁界Ｂｉ、地磁気Ｂｇ、及び外部磁界Ｂｘを測定する場合、複数の磁気データｑ１〜ｑＮの各々の表す座標は、内部磁界を表すベクトルＳＢｉ、地磁気を表すベクトルＳＢｇ（μ）、及び外部磁界を表すベクトルＳＢｘ（μ、Ｐｓ）の和を表すベクトルにより示される。従って、複数の磁気データｑ１〜ｑＮの示す座標は、中心点ｃＯＧを起点とする地磁気を表すベクトルＳＢｇ（μ）の終点を表す球面ＳＧと、中心点ｃＯＧを起点とする外部磁界を表すベクトルＳＢｘ（μ、Ｐｓ）の終点を表す曲面ＳＸとを、中心点ｃＯＧを起点として重ね合わせた立体ＳＤの表面近傍に分布する。
【０１１６】
外部磁界Ｂｘを表す曲面ＳＸが、球面とは異なる歪んだ形状を有する場合には、立体ＳＤも球面とは異なる歪んだ形状を有する。立体ＳＤが、球面とは異なる歪んだ形状を有する場合、複数の磁気データｑ１〜ｑＮの示す座標に基づいて、地磁気Ｂｇを表す球面ＳＧの中心点ｃＯＧの座標を算出することは困難である。仮に、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有するような球面Ｓを定め、球面Ｓの中心点ｃＳを算出しても、球面Ｓの中心点ｃＳと球面ＳＧの中心点ｃＯＧとは異なる座標となる可能性が高いからである。（図２４参照）。このように、立体ＳＤが球面とは異なる歪んだ形状を有し、球面ＳＧの中心点ｃＯＧの示す座標の算出が困難である場合には、複数の磁気データｑ１〜ｑＮに基づいたオフセットｃＯＦＦの算出を防止する必要がある。
【０１１７】
但し、不均一な外部磁界Ｂｘの影響が小さく、立体ＳＤの形状がほぼ球面と看做せる場合には、複数の磁気データｑ１〜ｑＮの示す座標に基づいて、球面ＳＧの中心点ｃＯＧの示す座標を算出することができる。例えば、図１８（Ａ）に示すように、外部磁界Ｂｘが微弱である場合、地磁気Ｂｇを表す球面ＳＧと、外部磁界Ｂｘを表す曲面ＳＸとの重ね合わせである立体ＳＤは、球面ＳＧとほぼ同一の図形となる。従って、複数の磁気データｑ１〜ｑＮの示す座標は、球面ＳＧの近傍に分布すると看做すことができるため、複数の磁気データｑ１〜ｑＮから、球面ＳＧの中心点ｃＯＧの示す座標を算出することができる。
また、図１８（Ｂ）に示すように、不均一な外部磁界Ｂｘが大きい場合であっても、立体ＳＤの形状がほぼ球面と看做せる場合がある。例えば、不均一な外部磁界Ｂｘが存在する場合であっても、Ｎ個の磁気データｑ１〜ｑＮを取得する際に、機器１または機器１ａの利用者が機器１または機器１ａを手で握って３次元磁気センサ６０の位置ＰＳが変化するように振るのではなく、３次元磁気センサ６０の位置ＰＳを固定して姿勢μのみを変化させる場合には、外部磁界Ｂｘは、センサ座標系ΣＳにおいて、姿勢μに基づいてその方向のみを変化させる一定の大きさのベクトルＳＢｘ（μ）として表現される。この場合、外部磁界Ｂｘを表す曲面ＳＸの形状は、中心点ｃＯＧを中心とする球面となるため、中心点ｃＯＧを中心とする球面ＳＧと、中心点ｃＯＧを中心とする球面の形状の曲面ＳＸとを、中心点ｃＯＧを中心に重ね合わせた立体ＳＤの形状は、中心点ｃＯＧを中心とする球面となる。従って、複数の磁気データｑ１〜ｑＮの示す座標に基づいて、立体ＳＤの表す球面の中心点の座標を算出することにより、球面ＳＧの中心点ｃＯＧの示す座標を算出することができる。
【０１１８】
本実施形態では、複数の磁気データｑ１〜ｑＮの示す座標に基づいて、外部磁界Ｂｘによる影響の大きさ、すなわち、立体ＳＤの形状と球面の形状との相違の程度を評価する。これにより、複数の磁気データｑ１〜ｑＮの示す座標に基づいて、オフセットｃＯＦＦが算出可能であるかを判断し、外部磁界Ｂｘの影響を受けた不正確なオフセットが算出されることを防止することが可能となる。
【０１１９】
なお、詳細は後述するが、本実施形態に係る地磁気測定装置は、楕円面補正部２００、楕円面球面変換部５００、及び歪判定部９００により、複数の磁気データｑ１〜ｑＮの示す座標の分布の形状と、楕円面の形状との相違の程度について評価することもできる。複数の磁気データｑ１〜ｑＮの示す座標を、楕円補正により複数の変換後磁気データｓ１〜ｓＮの表す座標に変換した場合、複数の変換後磁気データｓ１〜ｓＮの示す座標を近傍に有する立体ＳＤＥの形状が球面と看做すことができれば、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する立体ＳＤの形は、楕円面と看做すことができるからである。
すなわち、図１５に示すように、楕円面補正部２００、楕円面球面変換部５００、及び歪判定部９００は、複数の磁気データｑ１〜ｑＮの示す座標の分布の形状が、球面と看做すことのできる形状、楕円面と看做すことのできる形状、または、球面と看做すことも楕円面と看做すこともできない歪んだ形状、のいずれに該当するかを判断する、歪形状判断部４として機能する。
そして、歪形状判断部４が、複数の磁気データｑ１〜ｑＮの示す座標の分布の形状を、球面と看做すことのできる形状、または、楕円面と看做すことのできる形状であると判断した場合、地磁気測定装置は、オフセットｃＯＦＦを算出するが、球面とも楕円面とも異なる歪んだ形状であると判断した場合、オフセットｃＯＦＦを算出しない。
このように、本実施形態に係る地磁気測定装置は、歪形状判断部４を備えることにより、外部磁界Ｂｘの影響による不正確なオフセットが算出されることを防止することが可能になるとともに、外部磁界Ｂｘの影響が無視できる程度に小さい場合には、ソフトアイアン効果が生じている場合及び生じていない場合の双方において正確なオフセットｃＯＦＦの算出が可能になる。
【０１２０】
以下、本実施形態における、オフセットｃＯＦＦの候補となる座標の算出方法と、当該座標のオフセットｃＯＦＦへの採否の判断方法について、具体的に説明する。
【０１２１】
［８．  第２実施形態に係る地磁気測定装置のオフセット導出処理フロー］    
図１９は、第２実施形態に係る地磁気測定装置が、オフセットを導出する動作を説明するためのフローチャートである。このフローチャートは、ＣＰＵ１０が、本実施形態に係る磁気データ処理プログラムを実行することにより実施される。
【０１２２】
ステップＳ１において、地磁気測定装置は、初期化処理を行う。初期化処理は、蓄積部１００に備えられるバッファＢＵ１に格納された複数の磁気データｑ１〜ｑＮと、蓄積部１００に備えられるバッファＢＵ２に格納された各種データ（複数の変換後磁気データｓ１〜ｓＮ等）とを、破棄する処理である。なお、本実施形態の地磁気測定装置は、初期化処理において、バッファＢＵ１に格納された複数の磁気データｑ１〜ｑＮの全部を廃棄するが、古い方から一定割合のデータのみを廃棄してもよい。
【０１２３】
ステップＳ２において、地磁気測定装置は、磁気データ取得処理を行う。磁気データ取得処理は、３次元磁気センサ６０から順次出力される複数の磁気データｑ１〜ｑＮを、蓄積部１００に備えられるバッファＢＵ１に格納する処理である（Ｎは精度よくオフセットを導出するために必要な磁気データの規定測定回数を表す９以上の自然数）。
【０１２４】
ステップＳ３において、地磁気測定装置は、磁気データ分布判定処理を行う。磁気データ分布判定処理は、分布判定部７００が行う処理である。磁気データ分布判定処理において、分布判定部７００は、センサ座標系ΣＳにおいて、複数の磁気データｑ１〜ｑＮの示す座標の分布が、３次元的な広がりを有しているか否かを判定し、判定結果を出力する。
判定結果が肯定である場合、地磁気測定装置は、処理をステップＳ４に進める。一方、判定結果が否定である場合、すなわち、複数の磁気データｑ１〜ｑＮの示す座標の分布が２次元的または１次元的である場合、地磁気測定装置は、処理をステップＳ１に戻す。
【０１２５】
ステップＳ４において、地磁気測定装置は、中心点算出処理を行う。中心点算出処理は中心点算出部８００が行う処理である。中心点算出処理において、中心点算出部８００は、センサ座標系ΣＳにおいて、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する球面Ｓの中心点ｃＳの示す座標を算出し、これを出力する。
【０１２６】
ステップＳ５において、地磁気測定装置は、歪判定処理を行う。ステップＳ５における歪判定処理は、歪判定部９００が行う処理である。ステップＳ５の歪判定処理において、歪判定部９００は、複数の入力座標として複数の磁気データｑ１〜ｑＮの示す座標を適用したうえで、複数の入力座標がある立体ＳＤの近傍に分布すると仮定して、立体ＳＤの形状と球面の形状との相違の程度を評価することで、立体ＳＤの形状が球面と看做すことができるか否かを判定し、判定結果を出力する。
判定結果が肯定である場合、地磁気測定装置は、処理をステップＳ１０に進める。一方、判定結果が否定である場合、地磁気測定装置は、処理をステップＳ６に進める。
【０１２７】
ステップＳ６において、地磁気測定装置は、初期楕円面生成処理を行う。初期楕円面生成処理は、３節において説明した楕円面初期補正値生成部３００が行う処理である。前述の通り、初期楕円面生成処理において、楕円面初期補正値生成部３００は、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する初期楕円面ＶＥ０の中心点である初期中心点ｃＥ０の座標と、初期楕円面ＶＥ０上の座標を球面ＳＥ０上の座標へと変換する初期楕円面補正行列Ｔ０とを算出する。また、楕円面初期補正値生成部３００は、複数の磁気データｑ１〜ｑＮに基づいて算出される、第１楕円面係数行列Ｄｘｘ、第２楕円面係数行列Ｄｙｙ、第３楕円面係数行列Ｄｚｚ、中心点ｃｘｘの座標、中心点ｃｙｙの座標、及び中心点ｃｚｚの座標が、第１の条件及び第２の条件の双方を満たすか否かの判定を行う。当該判定結果が肯定である場合、地磁気測定装置は、処理をステップＳ７に進める。一方、判定結果が否定である場合、地磁気測定装置は、処理をステップＳ１に戻す。
なお、前述のとおり、楕円面初期補正値生成部３００は、第１の条件についての判定を行わず、第２の条件についてのみ判定を行うものであってもよい。
【０１２８】
ステップＳ７において、地磁気測定装置は最適楕円面生成処理を行う。最適楕円面生成処理は、４節において説明した楕円面最適補正値生成部４００の行う処理である。前述の通り、最適楕円面生成処理において、楕円面最適補正値生成部４００は、初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０に基づいて最適楕円面補正行列ＴＯＰ及び最適中心点ｃＥＯＰの座標を算出する。
【０１２９】
ステップＳ８において、地磁気測定装置は楕円面球面変換処理を行う。楕円面球面変換処理は、楕円面球面変換部５００が行う処理である。楕円面球面変換処理において、楕円面球面変換部５００は、最適楕円面補正行列ＴＯＰ及び最適中心点ｃＥＯＰに基づいて、最適楕円面ＶＥＯＰ近傍に存在する複数の磁気データｑ１〜ｑＮの示す座標を、複数の変換後磁気データｓ１〜ｓＮにより表される球面ＳＥＯＰ近傍の座標へと変換する。その後、楕円面球面変換部５００は、複数の変換後磁気データｓ１〜ｓＮを、蓄積部１００に備えられたバッファＢＵ２に格納する。バッファＢＵ２は、ＲＡＭ２０により構成される。
【０１３０】
ステップＳ９において、地磁気測定装置は、歪判定処理を行う。ステップＳ９における歪判定処理は、ステップＳ５における歪判定処理と同様に、歪判定部９００が行う処理である。ステップＳ９の歪判定処理において、歪判定部９００は、複数の入力座標として複数の変換後磁気データｓ１〜ｓＮの示す座標を適用したうえで、複数の入力座標が立体ＳＤＥの近傍に分布すると仮定して、立体ＳＤＥの形状と球面の形状との相違の程度を評価することで、立体ＳＤＥの形状が球面と看做すことができるか否かを判定し、判定結果を出力する。
判定結果が肯定である場合、地磁気測定装置は、処理をステップＳ１０に進める。一方、判定結果が否定である場合、地磁気測定装置は、処理をステップＳ１に戻す。
なお、以下では、ステップＳ５において行われる歪判定処理と、ステップＳ９において行われる歪判定処理とを区別する場合、前者を第１の歪判定処理と称し、後者を第２の歪判定処理と称する。また、第１の歪判定処理において形状が評価される立体ＳＤと区別するため、第２の歪判定処理において形状が評価される立体を立体ＳＤＥと称する。第１の歪判定処理と第２の歪判定処理とは、複数の入力座標の値が異なる以外は、同一の処理である。
【０１３１】
ステップＳ１０において、地磁気測定装置は、オフセット採用処理を行う。
オフセット採用処理は、オフセット採用部６１０ａが行う処理である。オフセット採用処理において、オフセット採用部６１０ａは、中心点ｃＳまたは中心点ｃＥＯＰの示す座標をオフセットｃＯＦＦとして採用するとともに、単位行列Ｉまたは最適楕円面補正行列ＴＯＰを楕円面補正行列ＴＥとして採用する。
具体的には、ステップＳ５における判定結果が肯定である場合、オフセット採用部６１０ａは、ステップＳ４において中心点算出部８００が算出した球面Ｓの中心点ｃＳの座標を表すベクトルをオフセットｃＯＦＦとして採用するとともに、単位行列Ｉを楕円面補正行列ＴＥとして採用する。一方、ステップＳ５における判定結果が否定であり、且つ、ステップＳ９における判定結果が肯定である場合、オフセット採用部６１０ａは、ステップＳ７で楕円面最適補正値生成部４００が算出した最適中心点ｃＥＯＰの座標を表すベクトルをオフセットｃＯＦＦとして採用するとともに、楕円面最適補正値生成部４００が算出した最適楕円面補正行列ＴＯＰを楕円面補正行列ＴＥとして採用する。そして、オフセット採用部６１０ａは、オフセットｃＯＦＦと楕円面補正行列ＴＥとを出力する。
また、ステップＳ９における判定結果が否定である場合、オフセット採用部６１０ａは、オフセットｃＯＦＦ及び楕円面補正行列ＴＥを採用しない。
なお、５節において説明したとおり、地磁気ベクトル計算部６２０は、オフセットｃＯＦＦと楕円面補正行列ＴＥとを用いて、３次元磁気センサ６０から出力される磁気データｑｉの示す座標に対して楕円補正を施すことで、地磁気Ｂｇの向きを算出する。地磁気ベクトル計算部６２０が楕円補正に用いるオフセットｃＯＦＦ及び楕円面補正行列ＴＥは、オフセット採用部６１０ａから出力されるオフセットｃＯＦＦ及び楕円面補正行列ＴＥにより更新される。
【０１３２】
本実施形態では、ステップＳ９において、判定結果が否定である場合、地磁気測定装置は、処理をステップＳ１に戻すが、その際に、何らかのメッセージを表示部５０に出力した上で処理を一旦停止させ、ユーザからの指示を待ってステップＳ１の処理を再開させてもよい。
Ｎ個の磁気データｑ１〜ｑＮを取得する際に、利用者が機器１を手で握って回転させるのではなく、機器１の位置を固定したままその姿勢のみを変化させるようにすると、外部磁界Ｂｘの影響を低く抑えることができる（図１８（Ｂ）参照）。そのため、ステップＳ９における判定結果が否定である場合には、機器１の位置を固定したまま回転させることを、利用者に対して指示してもよい。利用者に対する指示は、機器１の表示部５０に画像や動画等を表示したり、音声等を用いたりすることで、行ってもよい。
また、本実施形態では、ステップＳ６またはＳ９において、判定結果が否定である場合に、処理をステップＳ１に戻さずに、フローチャートに示された処理を終了してもよい。
【０１３３】
このように、本実施形態の地磁気測定装置は、ステップＳ５において、複数の磁気データｑ１〜ｑＮの示す座標が、球面Ｓの近傍に分布するか、球面とは異なる歪んだ形状の立体ＳＤの近傍に分布するかを判定する。そして、複数の磁気データｑ１〜ｑＮの示す座標が、球面とは異なる歪んだ形状の立体ＳＤの近傍に分布すると判定された場合、ステップＳ６〜Ｓ９を実行することにより、複数の磁気データｑ１〜ｑＮの示す座標が、最適楕円面ＶＥＯＰの近傍に存在するか否かの判定を行う。
すなわち、本実施形態の地磁気測定装置は、複数の磁気データｑ１〜ｑＮの示す座標の分布形状が、球面、楕円面、及び球面とも楕円面とも異なる歪んだ形状の立体のうち、いずれの形状と看做すことが適当であるかについて、判別することが可能となる。
これにより、本実施形態に係る地磁気測定装置は、複数の磁気データｑ１〜ｑＮの示す座標が、球面または楕円面の近傍に分布する場合には、これらの中心点をオフセットとして採用することにより、正確な地磁気の向きを算出することができる。一方、複数の磁気データｑ１〜ｑＮの示す座標が、球面とも楕円面とも異なる歪んだ形状の立体の近傍に分布する場合、本実施形態に係る地磁気測定装置は、オフセットの算出を防止するため、不正確な地磁気の向きの算出を防止することができる。
【０１３４】
以下において、磁気データ分布判定処理、中心点算出処理、及び歪判定処理の詳細を説明する。なお、理解を容易なものとするために、中心点算出処理を説明した後に、磁気データ分布判定処理の説明を行う。
【０１３５】
［９．  中心点算出処理］              
図２０を参照しつつ、ステップＳ４において、中心点算出部８００が行う中心点算出処理について説明する。中心点算出処理は、３次元磁気センサ６０が出力するＮ個の磁気データｑ１〜ｑＮの示す座標が、半径ｒＳの球面Ｓの近傍に分布すると仮定して、球面Ｓの中心点ｃＳの座標を算出する処理である。球面Ｓは、センサ座標系ΣＳ上で、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有するように定められる球面の中心点の座標を求めるために、計算の便宜上導入される球面であり、地磁気Ｂｇを表す球面ＳＧとは異なる球面である。なお、以下において記載されるベクトル及び座標は、特に断りがない場合には、センサ座標系ΣＳにおいて表現されたものであるとする。
【０１３６】
複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する球面Ｓの中心点ｃＳの座標の算出は、公知の方法を適宜適用すればよい。例えば、以下の方法により算出してもよい。
【０１３７】
磁気データｑｉの示す座標を式（１１）で表し、中心点ｃＳの座標を以下の式（７２）で表した場合、複数の磁気データｑ１〜ｑＮの示す座標が、半径ｒＳの球面Ｓ上に存在することは、以下の式（７１）で表される。
【数２５】
0005866893.tif 000026
【０１３８】
図２０に示すように、磁気データｑｉの示す座標を、複数の磁気データｑ１〜ｑＮの重心ｑＣの示す座標を起点とするベクトル（ｑｉ−ｑＣ）で表現した場合、式（７１）に対して複数の磁気データｑ１〜ｑＮの示す座標の各々を代入した式と、式（７１）とに基づいて、以下の式（７３）を得ることができる。以下では、式（７３）に示す方程式を、球面方程式と呼ぶ。ここで、重心ｑＣは、以下の式（７４）及び式（７５）で定義される３次元のベクトルである。また、行列Ｘは、式（７６）に示すＮ行３列の行列であり、ベクトルｊは、式（７７）に示すＮ次元のベクトルであり、値ＲＡＶＥは、式（７８）に示す値である。
【数２６】
0005866893.tif 000027
【０１３９】
式（７３）に示した球面方程式は、複数の磁気データｑ１〜ｑＮの示す座標の全てが、中心点ｃＳを中心とする球面Ｓ上に完全に一致する場合には解を有する。しかし、３次元磁気センサ６０の測定誤差等を考慮すると、複数の磁気データｑ１〜ｑＮの全てが、球面Ｓと完全に一致することはないため、球面方程式は解を持たない。そこで、統計的な手法により、球面方程式の尤もらしい解を得るために、式（７９）で表す誤差を吸収するベクトルである第１球面誤差ベクトルδＳを導入する。ここで、式（７９）に現れる変数ベクトルｃは、式（５）に示した３次元のベクトルであるが、本節では、中心点ｃＳの座標を表すための変数として用いる。
【数２７】
0005866893.tif 000028
【０１４０】
第１球面誤差ベクトルδＳのノルムを最小にするベクトルｃ、換言すれば、(δＳ)Ｔ(δＳ)を最小化するようなベクトルｃの表す座標が、球面Ｓの中心点ｃＳの座標として尤もらしいものであるといえる。ここで、以下の式（８０）で示される中心点算出関数ｆＳ（ｃ）を定義すると、中心点算出関数ｆＳ（ｃ）を最小化するベクトルｃの表す座標が、球面Ｓの中心点ｃＳの座標として尤もらしい値となる。中心点ｃＳの座標は、式（８２）に示す３行３列の共分散行列Ａが正則である場合には、式（８１）により算出される。
【数２８】
0005866893.tif 000029
【０１４１】
前述の通り、３次元磁気センサ６０が、内部磁界Ｂｉ及び地磁気Ｂｇのみを検出する場合には、球面Ｓと、地磁気Ｂｇを表す球面ＳＧとは、ほぼ同一の球面となり、球面Ｓの中心点ｃＳと、球面ＳＧの中心点ｃＯＧとは、ほぼ同一の座標となる。従って、３次元磁気センサ６０が内部磁界Ｂｉ及び地磁気Ｂｇのみを検出する場合、式（８１）に示される中心点ｃＳの座標を示すベクトルを、磁気センサのオフセットｃＯＦＦとして採用することができる。
【０１４２】
［１０．  磁気データ分布判定処理］                                        
以下では、ステップＳ３において、分布判定部７００が行う磁気データ分布判定処理について説明する。
【０１４３】
前述した、中心点算出処理において、球面Ｓの中心点ｃＳを算出するためには、複数の磁気データｑ１〜ｑＮの示す座標が、センサ座標系ΣＳにおいて３次元的な広がりを有して分布していることを必要とする。しかし、機器１（３次元磁気センサ６０）の姿勢μは、機器１の利用者が手で機器１を動かすことにより変化するため、機器１の動かし方が不十分な場合には、機器１の姿勢変化は３次元的とはならず２次元的なものとなることがある。この場合、センサ座標系ΣＳにおける複数の磁気データｑ１〜ｑＮの示す座標は、３次元的な広がりを有さず２次元的に分布する。
例えば、図２１に示すように、複数の磁気データｑ１〜ｑＮの示す座標が、センサ座標系ΣＳの平面π上の円πＣ近傍に２次元的に分布する場合、球面Ｓは、円πＣを切断面に有するような球面であるということしか特定できない。円πＣを切断面に有するような球面は、円πＣの中心点πＣＯを通り平面πに直交する直線πＬ上の中心点ｃπ１を中心とする球面Ｓπ１であるかもしれないし、直線πＬ上の中心点ｃπ２を中心とする球面Ｓπ２であるかもしれない。つまり、球面Ｓの中心点ｃＳは、直線πＬ上に位置することまでは特定可能であるが、具体的に直線πＬ上のどの位置に存在するかについては特定不可能である。このように、複数の磁気データｑ１〜ｑＮの示す座標が２次元的に分布する場合、複数の磁気データｑ１〜ｑＮに基づいて、正確な中心点ｃＳを算出することはできない。
【０１４４】
複数の磁気データｑ１〜ｑＮに基づいて球面Ｓの中心点ｃＳを算出するためには、図２２に示すように、センサ座標系ΣＳにおいて複数の磁気データｑ１〜ｑＮの示す座標が３次元的に広がりを持って分布していることが必要となる。磁気データ分布判定処理において、分布判定部７００は、複数の磁気データｑ１〜ｑＮの示す座標が３次元的に分布しているか否かを判定する。複数の磁気データｑ１〜ｑＮの示す座標が３次元的に分布しているか否かの判定は、公知の方法を適宜適用すればよいが、例えば、式（８２）に示される共分散行列Ａを用いて判定してもよい。以下において、共分散行列Ａの性質を説明する。
【０１４５】
共分散行列Ａの固有値を、大きい順番に最大固有値λ１、中間固有値λ２、最小固有値λ３とし、それぞれの固有値に対応する大きさ１に正規化された固有ベクトルをｕ１、ｕ２、ｕ３とする。また、磁気データｑｉを、前述した重心ｑＣを原点とする重心座標系ΣＣにおいて表したベクトルをＣｑｉと表す。このとき、固有値λｊ（ｊ＝１、２、３）は、固有ベクトルｕｊ方向の分散σ２ｊに等しい。
図２２に示すように、各固有ベクトルｕ１〜ｕ３を重心座標系ΣＣの原点ｑＣを起点となるように配置する。このとき、例えばｊ＝１の場合について検討する。固有値λ１は、ベクトルＣｑｉを、固有ベクトルｕ１へ射影した長さＬｉ１の二乗（Ｌｉ１）２を、Ｎ個の磁気データＣｑｉ（ｉ＝１、２、…Ｎ）について平均した値に等しい。すなわち、固有値λｊは、Ｎ個の磁気データＣｑｉが、重心ｑＣから固有ベクトルｕｊ方向にどの程度離れているか、つまり、複数の磁気データｑ１〜ｑＮの示す座標の分布が固有ベクトルｕｊの方向にどの程度の広がり有するかを表す。
【０１４６】
最小固有値λ３に対応する固有ベクトルｕ３の方向が、複数の磁気データｑ１〜ｑＮの示す座標の分布の広がりが最も小さい方向であり、最小固有値λ３が、複数の磁気データｑ１〜ｑＮの示す座標の分布の広がりが最も小さい方向における広がりの程度を示す指標である。従って、複数の磁気データｑ１〜ｑＮの示す座標が３次元的に分布しているといえるためには、最小固有値λ３の値が、一定の閾値（分散許容値）λ０以上であればよい。
【０１４７】
磁気データ分布判定処理において、分布判定部７００は、共分散行列Ａの最小固有値λ３が閾値λ０以上であれば、複数の磁気データｑ１〜ｑＮの示す座標の分布が十分に３次元的であると判断して、処理を前述したステップＳ４の中心点算出処理へと進める。一方、分布判定部７００は、最小固有値λ３が閾値λ０未満である場合には、複数の磁気データｑ１〜ｑＮの示す座標が３次元的な広がりを有さないと判断し、処理をステップＳ１の初期化処理に戻す。
【０１４８】
［１１．  歪判定処理］                                                  
歪判定部９００は、ステップＳ５において第１の歪判定処理を行うとともに、ステップＳ９において第２の歪判定処理を行う。第２の歪判定処理は、複数の入力座標として、複数の磁気データｑ１〜ｑＮの示す座標の代わりに、複数の変換後磁気データｓ１〜ｓＮの示す座標を用いて歪判定処理を行う点を除いて、第１の歪判定処理と同様の処理である。
以下１１．１節において、第１の歪判定処理について述べた上で、１１．２節において、第２の歪判定処理について述べる。
【０１４９】
［１１．１．  第１の歪判定処理］                                          
歪判定処理は、複数の入力座標、すなわち、複数の磁気データｑ１〜ｑＮの示す座標が、球面とは異なる歪んだ形状の立体ＳＤの表面近傍に分布すると仮定する。立体ＳＤは、図２３に示すように、球面（第２球面）Ｓ２と、歪誤差ベクトルｋ（Ｅ）とを、足し合わせることで得られる図形であり、以下の式（８３）により表される。以下では、式（８３）に示す方程式を、立体方程式と呼ぶ。
ここで、球面Ｓ２は、中心点（第２球面の中心点）ｃＳ２を中心とする球面であり、立体方程式のうち、歪誤差ベクトルｋ（Ｅ）以外の成分「Ｘ（ｃ−ｑＣ）−ｊ」として表される。
歪誤差ベクトルｋ（Ｅ）は、以下の式（８４）により表されるＮ次元ベクトルである。但し、歪評価行列Ｅは、以下の式（８５）に表される３行３列の対称行列であり、基準点ｗＫＥは、以下の式（８６）により表される３次元ベクトルである。また、式（８３）の右辺に現れる０ＮはＮ次元の零ベクトルである。式（８３）の左辺に現れる変数ベクトルｃは、式（５）で示された３次元のベクトルであるが、本節では、球面Ｓ２の中心点ｃＳ２を表すための変数として用いる。
歪判定処理は、立体方程式のうち、歪を表す成分ｋ（Ｅ）の大きさを評価することにより、立体ＳＤの形状と、球面Ｓ２の形状とが、どの程度相違しているかについて評価する。具体的には、立体方程式における歪誤差ベクトルｋ（Ｅ）の影響の大きさを、後述する式（９３）及び式（９４）で示される歪評価値ｇＤ（Ｅ）に基づいて評価する。
【数２９】
0005866893.tif 000030
【０１５０】
Ｎ次元の歪誤差ベクトルｋ（Ｅ）を構成するＮ個の要素のうち、第ｉ行目の要素ｋｅ(ｑｉ−ｗＫＥ)は、以下の式（８７）で表される関数ｋｅ（ｖ）に対して、基準点ｗＫＥの示す座標を起点として磁気データｑｉの示す座標を表すベクトル（ｑｉ−ｗＫＥ）を代入することにより与えられる。関数ｋｅ（ｖ）は、式（８５）に示す歪評価行列Ｅを係数行列とし、式（８８）に示すベクトルｖの３つの要素を変数とする、二次形式で表現される関数である。つまり、関数ｋｅ（ｖ）は、ベクトルｖと、ベクトルｖを歪評価行列Ｅにより変換したベクトルＥｖとの内積を示す。
なお、第１の歪判定処理では、基準点ｗＫＥとして、以下の式（８９）に示すように、球面Ｓの中心点ｃＳを採用する。
【数３０】
0005866893.tif 000031
【０１５１】
３次元磁気センサ６０の測定誤差等を考慮すると、複数の磁気データｑ１〜ｑＮの示す座標の全てが、立体ＳＤと完全に一致する位置に存在することはないため、式（８３）に示した立体方程式は解を持たない。そこで、統計的な手法により、立体方程式の解として尤もらしい値を得るために、式（９０）で表される誤差を吸収するベクトルである立体誤差ベクトルδＳＤを導入する。立体誤差ベクトルδＳＤは、以下の式（９１）に示す第２球面誤差ベクトルδＳ２と、歪誤差ベクトルｋ（Ｅ）とを加算したものである。第２球面誤差ベクトルδＳ２は、立体方程式のうち、球面Ｓ２を表す成分「Ｘ（ｃ−ｑＣ）−ｊ」に対応する成分である。
立体誤差ベクトルδＳＤは、複数の磁気データｑ１〜ｑＮの示す座標と、立体ＳＤの表面との誤差を表す、Ｎ次元のベクトルである。立体誤差ベクトルδＳＤのノルムを最小にする変数ベクトルｃ及び歪評価行列Ｅ、換言すれば、以下の式（９２）で示す歪評価関数ｆＳＤ（Ｅ、ｃ）を最小化する変数ベクトルｃ及び歪評価行列Ｅによって、複数の磁気データｑ１〜ｑＮの示す座標を表面近傍に有する立体ＳＤが表現される。
【数３１】
0005866893.tif 000032
【０１５２】
以下で、式（９０）に示した立体誤差ベクトルδＳＤの性質を、式（７９）に示した第１球面誤差ベクトルδＳの性質と対比しつつ説明する。
まず、第１球面誤差ベクトルδＳは、複数の磁気データｑ１〜ｑＮの示す座標と、球面Ｓとの誤差を吸収するためのベクトルである。第１球面誤差ベクトルδＳを構成する１行目〜Ｎ行目の各要素は、それぞれ独立した変数である。従って、第１球面誤差ベクトルδＳによって、複数の磁気データｑ１〜ｑＮが示す座標と球面Ｓとの誤差を吸収する場合には、複数の磁気データｑ１〜ｑＮが示す座標と球面ＳとのＮ個の誤差の各々は、互いに制約の無い、Ｎ個が互いに独立して定められる値となる。つまり、第１球面誤差ベクトルδＳによって表されるＮ個の誤差は、それぞれが独立に確率的に定められるものであり、Ｎ個の誤差は全体として、対称性を有し、かつ方向依存性のないホワイトノイズである。
すなわち、中心点算出処理は、複数の磁気データｑ１〜ｑＮの示す座標と、球面Ｓとの誤差をホワイトノイズである第１球面誤差ベクトルδＳにより表現し、第１球面誤差ベクトルδＳを最小化する球面Ｓの中心点ｃＳを求める処理である。
【０１５３】
一方、立体誤差ベクトルδＳＤは、第２球面誤差ベクトルδＳ２と歪誤差ベクトルｋ（Ｅ）との和によって表されるベクトルであり、磁気データｑ１〜ｑＮの示す座標と立体ＳＤとの誤差を吸収するベクトルである。
第２球面誤差ベクトルδＳ２は、第１球面誤差ベクトルδＳと同様に、磁気データｑ１〜ｑＮの示す座標と、球面Ｓ２との誤差を、ホワイトノイズとして表現するベクトルである。
一方、歪誤差ベクトルｋ（Ｅ）は、式（８７）で示した３変数二次形式の関数ｋｅ（ｖ）を各要素とするベクトルである。３変数の二次形式は、変数が２次の項から構成される関数であり、３次元空間上の様々な曲面、例えば、直線、平面、柱面、球面、楕円面、錐面、１葉双曲面、２葉双曲面、及び各種放物面等を描くことができる。従って、歪誤差ベクトルｋ（Ｅ）は、各磁気データｑ１〜ｑＮが示す座標と球面Ｓ２とのＮ個の誤差の各々を、互いに独立したものではなく、Ｎ個の誤差の全てが同一の関数ｋｅ（ｖ）により表される３次元空間上の曲面上に存在するという制約を持った値として表現する。
このように、立体誤差ベクトルδＳＤは、各磁気データｑ１〜ｑＮが示す座標と球面Ｓ２とのＮ個の誤差を、ホワイトノイズである第２球面誤差ベクトルδＳ２と、球面Ｓ２からの歪を表す曲面を表す歪誤差ベクトルｋ（Ｅ）とに分離して表現する。
【０１５４】
立体方程式における歪誤差ベクトルｋ（Ｅ）の影響が、無視できる程度に小さい場合、立体ＳＤと球面Ｓ２とは、等しい図形と看做すことができ、式（９２）に定める歪評価関数ｆＳＤ（Ｅ、ｃ）と、式（８０）に定める中心点算出関数ｆＳ（ｃ）とは、等しい関数と看做すことができる。このとき、歪評価関数ｆＳＤ（Ｅ、ｃ）を最小化して求められる立体ＳＤと、中心点算出関数ｆＳ（ｃ）を最小化して求められる球面Ｓとは、等しいものと看做すことができるため、立体ＳＤの表面近傍に分布する複数の磁気データｑ１〜ｑＮが示す座標は、球面Ｓの近傍にも分布すると看做すことができる。前述の通り、立体ＳＤの形状が球面と看做せる場合、立体ＳＤが表す球面の中心点の示す座標と、球面ＳＧの中心点ｃＯＧとは、一致すると看做すことができる。よって、球面Ｓの中心点ｃＳの示す座標と、球面ＳＧの中心点ｃＯＧの示す座標は、等しいと看做すことができる。
このように、立体方程式における歪誤差ベクトルｋ（Ｅ）の影響が小さい場合、中心点算出手段において算出した球面Ｓの中心点ｃＳの示す座標は、球面ＳＧの中心点ｃＯＧの示す座標と等しいと看做されるため、中心点ｃＳの示す座標をオフセットｃＯＦＦとして採用することができる。
【０１５５】
逆に、立体方程式における歪誤差ベクトルｋ（Ｅ）の影響が大きい場合、図２３に示すように、複数の磁気データｑ１〜ｑＮが示す座標と球面Ｓ２との誤差は、ホワイトノイズである第２球面誤差ベクトルδＳ２と、球面Ｓ２からの歪を表す歪誤差ベクトルｋ（Ｅ）とによって吸収される。この場合、立体ＳＤの形状は、球面とは異なる形状を有する。
また、立体方程式における歪誤差ベクトルｋ（Ｅ）の影響が大きい場合、歪評価関数ｆＳＤ（Ｅ、ｃ）と、中心点算出関数ｆＳ（ｃ）とは、異なる関数となる。この場合、図２４に示すように、歪評価関数ｆＳＤ（Ｅ、ｃ）を最小化して求められる立体ＳＤと、中心点算出関数ｆＳ（ｃ）を最小化して求められる球面Ｓとは、異なる図形となり、立体ＳＤの表面近傍に分布する複数の磁気データｑ１〜ｑＮの示す座標は、球面Ｓの近傍に分布すると看做すことはできない。
中心点算出処理は、複数の磁気データｑ１〜ｑＮの示す座標が球面Ｓの近傍に存在していることを前提として、球面ＳＧの中心点ｃＯＧと等しいと看做すことのできる中心点ｃＳの示す座標を算出する処理である。従って、複数の磁気データｑ１〜ｑＮの示す座標が球面Ｓの近傍に存在しない場合、中心点ｃＳと、中心点ｃＯＧとは、一致しない。この場合、中心点ｃＳの座標を示すベクトルを、オフセットｃＯＦＦとして採用することはできない。
【０１５６】
このように、立体方程式における歪誤差ベクトルｋ（Ｅ）の影響の大きさを評価することにより、球面Ｓの中心点ｃＳを、オフセットｃＯＦＦとして採用することの可否を判断することが可能である。以下では、立体方程式おける歪誤差ベクトルｋ（Ｅ）の影響の大きさを評価する方法について述べる。
【０１５７】
ここで、以下の式（９３）及び（９４）に示す歪評価値ｇＤ（Ｅ）を、立体方程式における歪誤差ベクトルｋ（Ｅ）の影響の大きさを評価する評価値として定義する。歪評価値ｇＤ（Ｅ）は、歪評価行列Ｅの有する３つの固有値のうち、絶対値が最大となる最大固有値λＥ１の絶対値（つまり、歪評価行列Ｅのノルム）である。
歪評価値ｇＤ（Ｅ）が、歪許容値δ０以下の小さな値であれば、立体ＳＤと球面Ｓ２とは等しい図形であると看做すことができ、立体ＳＤの表面近傍に分布する複数の磁気データｑ１〜ｑＮの示す座標は、球面Ｓの近傍に分布すると看做すこともできる。この場合、中心点算出処理により求められた球面Ｓの中心点ｃＳの座標を表すベクトルを磁気センサのオフセットｃＯＦＦとして採用することができる。
【数３２】
0005866893.tif 000033
【０１５８】
ところで、前述の通り、歪誤差ベクトルｋ（Ｅ）の各要素は、基準点ｗＫＥの示す座標から見た磁気データｑｉの示す座標を表すベクトル（ｑｉ−ｗＫＥ）と、当該ベクトル（ｑｉ−ｗＫＥ）を歪評価行列Ｅにより変換したベクトルＥ（ｑｉ−ｗＫＥ）との内積である。
つまり、歪誤差ベクトルｋ（Ｅ）を構成する複数の要素の絶対値は、当該要素に対応する磁気データｑｉの示す座標を基準点ｗＫＥの示す座標から表したベクトル（ｑｉ−ｗＫＥ）と、歪評価行列Ｅが有する３つの固有値のうち絶対値の大きさが最大となる最大固有値λＥ１に対応する固有ベクトルｕＥ１とが平行となる場合には、大きな値となる。
従って、球面Ｓ２との誤差が大きな座標を示す磁気データｑｉが多く存在する領域を、基準点ｗＫＥの示す座標から表した方向と、歪評価行列Ｅの最大固有値λＥ１に対応する固有ベクトルｕＥ１の方向とが等しくなるように歪評価行列Ｅの各成分を定める場合、歪誤差ベクトルｋ（Ｅ）は、複数の磁気データｑ１〜ｑＮの示す座標と、球面Ｓ２との誤差の大きさを正確に表現する。
歪評価関数ｆＳＤ（Ｅ、ｃ）を最小化する歪評価行列Ｅは、複数の磁気データｑ１〜ｑＮの示す座標と、球面Ｓ２との誤差を正確に表現するように定められる。従って、歪評価行列Ｅの各成分は、最大固有値λＥ１に対応する固有ベクトルｕＥ１の向きと、球面Ｓ２からの誤差が大きな磁気データが多く存在する領域を基準点ｗＫＥの示す座標から表したベクトルの方向とが、近づくように定められる。そして、歪評価行列Ｅの最大固有値λＥ１は、球面Ｓ２からの誤差が大きな磁気データｑｉについての、その誤差の大きさを表現する値となる。
本実施形態は、立体ＳＤと球面との形状の相違の程度を示す歪評価値ｇＤ（Ｅ）を、歪評価行列Ｅの最大固有値λＥ１の値に基づいて定める。従って、歪評価値ｇＤ（Ｅ）によって、球面Ｓ２からの距離が大きな座標を示す磁気データｑｉと、球面Ｓ２との誤差の大きさ、すなわち、立体ＳＤと球面との形状の相違の程度を評価することが可能となる。
【０１５９】
以下に、歪評価値ｇＤ（Ｅ）を求める方法について説明する。
まず、式（８７）に示した関数ｋｅ（ｖ）は、以下の式（９５）に変形できる。また、歪誤差ベクトルｋ（Ｅ）の第ｉ行目の要素ｋｅ(ｑｉ−ｗＫＥ)は、式（９７）によって示される６次元のベクトルｋｅ２（ｉ）、及び、式（９８）によって表される歪評価行列Ｅの各成分を並べた６次元のベクトルｅＥを用いて、以下の式（９６）に変形できる。
【数３３】
0005866893.tif 000034
【０１６０】
ここで、式（９９）で示す行列Ｘ２を導入する。行列Ｘ２は、ベクトルｋｅ２（ｉ）を転置した１行６列のベクトルと、ベクトル（ｑｉ−ｑＣ）を転置した１行３列のベクトルとを、各行に並べて生成された、Ｎ行９列の行列である。
【数３４】
0005866893.tif 000035
【０１６１】
行列Ｘ２を用いて、式（９２）に示した歪評価関数ｆＳＤ（Ｅ、ｃ）は、以下の式（１００）に表す関数ｇＳＤ（ｅ）に変形される。なお、ベクトルｅは、以下の式（１０１）に示すように、ベクトルｅＥと、以下の式（１０２）に示す３次元のベクトルｅＸとを並べた、９次元のベクトルである。
【数３５】
0005866893.tif 000036
【０１６２】
式（１００）に示す関数ｇＳＤ（ｅ）を最小化する解ｅ＝ｅ０は、以下の式（１０３）に示す連立方程式に対して、ガウス消去法や、コレスキー分解法などを適用することで求められる。なお、式（１０３）は、式（１００）に対して、最小二乗法を適用することによって算出される正規方程式である。
【数３６】
0005866893.tif 000037
【０１６３】
このようにして得た解ｅ０に基づき、式（８５）の歪評価行列Ｅを復元する。そして、式（９３）に示した歪評価値ｇＤ（Ｅ）の値、即ち、歪評価行列Ｅのノルムを求め、歪評価値ｇＤ（Ｅ）の値が歪許容値δ０以下であるか否かを判定する。なお、歪評価行列Ｅのノルムは、歪評価行列Ｅの有する３つの固有値のうち絶対値が最大となる固有値λＥ１の絶対値に等しいため、ヤコビ法や、冪乗法により求めることができる。
【０１６４】
歪評価値ｇＤ（Ｅ）が歪許容値δ０以下である場合、地磁気測定装置は、処理をステップＳ１０のオフセット採用処理に進め、球面Ｓの中心点ｃＳの座標を表すベクトルを、オフセットｃＯＦＦとして採用する。
一方、歪評価値ｇＤ（Ｅ）が歪許容値δ０よりも大きな値である場合には、球面Ｓの中心点ｃＳの座標を表すベクトルを、オフセットｃＯＦＦとして採用することはできない。この場合、地磁気測定装置は、処理をステップＳ６の初期楕円面生成処理に進める。
【０１６５】
このように、第１の歪判定処理は、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する立体ＳＤと球面との形状の相違の程度を評価する。立体ＳＤと球面との形状の相違が無視できる程度に小さい場合には、立体ＳＤを球面と看做すことができるため、球面Ｓの中心点ｃＳの座標を表すベクトルを、オフセットｃＯＦＦとして採用することができる。この場合には、ソフトアイアン効果が生じていないと看做すことができるため、ステップＳ６〜Ｓ８に示す楕円補正（初期楕円面生成処理、最適楕円面生成処理、及び楕円面球面変換処理）を行うことなく、地磁気Ｂｇの向きを算出することができる。
すなわち、本実施形態に係る地磁気測定装置は、第１の歪判定処理を行うことにより、ソフトアイアン効果の有無を判断し、ソフトアイアン効果が生じていないと判断される場合には、楕円補正を行うことなく地磁気Ｂｇの向きを算出する。これにより、本実施形態に係る地磁気測定装置は、地磁気Ｂｇの向きの算出に係る計算負荷を大幅に低減することが可能となる。
【０１６６】
なお、本実施形態では、歪評価値ｇＤ（Ｅ）が歪許容値δ０以下である場合に、オフセット採用部６１０ａは、球面Ｓの中心点ｃＳの座標を表すベクトルをオフセットｃＯＦＦとして採用するが、球面Ｓ２の中心点ｃＳ２の座標を表すベクトルをオフセットｃＯＦＦとして採用してもよい。歪評価値ｇＤ（Ｅ）の値が閾値δ０以下である場合には、球面Ｓの中心点ｃＳの示す座標と、球面Ｓ２の中心点ｃＳ２の示す座標とは、ほぼ等しくなり、いずれもオフセットｃＯＦＦとして採用することができるからである。
なお、球面Ｓ２の中心点ｃＳ２の座標は、関数ｇＳＤ（ｅ）を最小化する解ｅ０のうち式（１０１）のｅＸに相当する３次元のベクトルを、式（１０２）に代入した場合の、変数ベクトルｃとして算出される。
【０１６７】
［１１．２．  第２の歪判定処理］  
図２５を参照しつつ、ステップＳ９において、歪判定部９００が行う第２の歪判定処理について説明する。
第２の歪判定処理は、ステップＳ５における第１の歪判定処理における判定結果が否定である場合、すなわち、図２５（Ａ）に示すように、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する立体ＳＤの形状が球面とは異なる歪んだ形状であると判定された場合に、図２５（Ｂ）に示すように、複数の変換後磁気データｓ１〜ｓＮの示す座標の分布形状を評価するための処理である。すなわち、第２の歪判定処理において、歪判定部９００は、複数の入力座標として、複数の変換後磁気データｓ１〜ｓＮの示す座標を用いる。
【０１６８】
第１の歪判定処理において、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する立体ＳＤの形状が球面とは異なる歪んだ形状であると評価される場合、立体ＳＤの歪の原因が、不均一な外部磁界Ｂｘではなく、ソフトアイアン効果である場合がある。不均一な外部磁界Ｂｘは存在せず、且つ、ソフトアイアン効果が生じている場合、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する立体ＳＤは、楕円面ＶＥと同一の形状を有すると看做すことができる。この場合、複数の磁気データｑ１〜ｑＮの示す座標を、最適楕円面補正行列ＴＯＰにより変換して得られる複数の変換後磁気データｓ１〜ｓＮの示す座標は、球面ＳＥＯＰ近傍に分布するため、複数の変換後磁気データｓ１〜ｓＮの示す座標を近傍に有する立体ＳＤＥの形状は、球面と看做すことができる。
逆に、不均一な外部磁界Ｂｘが存在する場合、図２５（Ａ）に示すように、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する立体ＳＤは、球面とは異なる歪んだ形状を有し、且つ、楕円面とも異なる形状を有する。この場合、複数の磁気データｑ１〜ｑＮの示す座標を、最適楕円面補正行列ＴＯＰにより変換して得られる複数の変換後磁気データｓ１〜ｓＮの示す座標は、図２５（Ｂ）に示すように球面ＳＥＯＰとは異なる歪んだ形状の立体ＳＤＥの近傍に分布する。
【０１６９】
このように、第２の歪判定処理は、複数の変換後磁気データｓ１〜ｓＮの示す座標を近傍に有する立体ＳＤＥの形状と球面（例えば、球面ＳＥＯＰ）の形状との相違の程度を評価することにより、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する立体ＳＤの形状と楕円面（例えば、楕円面ＶＥ）の形状との相違の程度を評価する処理である。
【０１７０】
なお、第２の歪判定処理において、立体ＳＤＥの形状が、球面の形状と等しいと看做すことができる場合、ソフトアイアン効果が生じているが、不均一な外部磁界Ｂｘの影響は存在しないため、最適中心点ｃＥＯＰの示す座標をオフセットｃＯＦＦとして採用することができる。
一方、図２５（Ｂ）に示すように、第２の歪判定処理において、立体ＳＤＥの形状が、球面の形状とは異なる歪を持った形状であると評価される場合、複数の磁気データｑ１〜ｑＮは不均一な外部磁界Ｂｘの影響を受けたものであるため、複数の磁気データｑ１〜ｑＮに基づいてオフセットｃＯＦＦの算出をすることはできない。
【０１７１】
以下において、第２の歪判定処理の具体的な処理内容を説明する。
前述の通り、第２の歪判定処理は、歪判定処理に用いる複数の入力座標として、複数の磁気データｑ１〜ｑＮの示す座標を用いる代わりに、複数の変換後磁気データｓ１〜ｓＮの示す座標を用いる点を除いて、第１の歪判定処理と同様の処理である。つまり、第２の歪判定処理は、第１の歪判定処理で用いた複数の磁気データｑ１〜ｑＮの示す座標に対して、以下の式（１０４）、式（１０５）に示す複数の変換後磁気データｓ１〜ｓＮの示す座標の値を代入した後に、１１．１節で述べた第１の歪判定処理を実行する処理である。
なお、歪判定処理で用いる値のうち、９節で述べた中心点算出処理において算出される値、例えば、式（８３）に示す立体方程式に現れる行列Ｘ等も、複数の磁気データｑ１〜ｑＮの示す座標代わりに、複数の変換後磁気データｓ１〜ｓＮの示す座標を用いて算出される。例えば、式（７４）に示す重心ｑＣには、以下の式（１０６）に示す、複数の変換後磁気データｓ１〜ｓＮの示す座標の重心ｓＣが代入されたうえで、第２の歪判定処理が行われる。これらの値の算出は中心点算出部８００において行われるが、歪判定部９００が行ってもよい。
また、基準点ｗＫＥには、中心点ｃＳの示す座標の代わりに、以下の式（１０７）に示す最適中心点ｃＥＯＰの示す座標を代入して、第２の歪判定処理が行われる。
第２の歪判定処理は、このようにして複数の変換後磁気データｓ１〜ｓＮの示す座標に基づいて定められる歪評価関数ｆＳＤ（Ｅ、ｃ）の値を最小化して求められる歪評価行列Ｅから歪評価値ｇＤ（Ｅ）を算出することにより、立体ＳＤＥと球面との形状の相違の程度を評価する。
【数３７】
0005866893.tif 000038
【０１７２】
複数の変換後磁気データｓ１〜ｓＮの示す座標に基づき算出された歪評価値ｇＤ（Ｅ）が歪許容値δ０以下である場合、ソフトアイアン効果が生じているが、不均一な外部磁界Ｂｘの影響は存在しないため、地磁気測定装置は、処理をステップＳ１０のオフセット採用処理に進め、最適中心点ｃＥＯＰの座標を表すベクトルを、オフセットｃＯＦＦとして採用する。前述の通り、第２の歪判定処理における球面Ｓ２の中心点ｃＳ２の座標を表すベクトルをオフセットｃＯＦＦとして採用してもよい。
一方、複数の変換後磁気データｓ１〜ｓＮの示す座標に基づき算出された歪評価値ｇＤ（Ｅ）が歪許容値δ０よりも大きな値である場合、不均一な外部磁界Ｂｘの影響が存在するため、地磁気測定装置は、処理をステップＳ１の初期化処理に戻し、球面ＳＥＯＰの最適中心点ｃＥＯＰの座標を表すベクトルが、オフセットｃＯＦＦとして採用されることを防止する。
なお、本実施形態において、第１の歪判定処理における歪許容値δ０と、第２の歪判定処理における歪許容値δ０とは等しい値に設定されるが、異なる値に設定してもよい。
【０１７３】
［１２．  第２実施形態の結論］                    
以上に示したように、第２実施形態に係る地磁気測定装置は、歪判定部９００を備えることにより、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する立体ＳＤと球面との形状の相違の程度を評価する。
そして、立体ＳＤの形状が球面と看做すことができる場合、簡易な計算により地磁気Ｂｇの向きを算出することができる。具体的には、第２実施形態に係る地磁気測定装置は、複数の磁気データｑ１〜ｑＮの示す座標が球面近傍に分布すると歪判定部９００が判定した場合、中心点算出部８００が算出する中心点ｃＳの座標を表すベクトルを、オフセットｃＯＦＦとして採用する。そして、地磁気測定装置は、中心点ｃＳの示す座標と、磁気データｑｉの示す座標とに基づいて、地磁気Ｂｇの向きを算出する。
このように、第２実施形態に係る地磁気測定装置は、３次元磁気センサ６０が軟磁性材料を備えない機器１ａに搭載され、ソフトアイアン効果が生じていない場合には、楕円補正を行うこと無く地磁気Ｂｇの向きを算出するため、計算負荷を軽減することができる。
【０１７４】
また、第２実施形態に係る地磁気測定装置は、楕円面補正部２００、楕円面球面変換部５００、及び歪判定部９００を備える。そして、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する立体ＳＤの形状が球面とは異なる歪んだ形状であると歪判定部９００が判定した場合に、歪判定部９００は、楕円面球面変換部５００の算出する複数の変換後磁気データｓ１〜ｓＮの示す座標を近傍に有する立体ＳＤＥと球面との形状の相違の程度を評価することで、立体ＳＤの歪の原因が、ソフトアイアン効果によるものであるのか、または、不均一な外部磁界Ｂｘによるものであるのか、を区別する。
そして、立体ＳＤＥの形状が球面とは異なる歪んだ形状である場合、つまり、立体ＳＤの歪の原因が不均一な外部磁界Ｂｘである場合には、地磁気測定装置は、不均一な外部磁界Ｂｘの影響下で測定された複数の磁気データｑ１〜ｑＮの示す座標に基づいてオフセットｃＯＦＦが算出されることを防止する。
一方、立体ＳＤＥの形状が球面であると看做せる場合、つまり、ソフトアイアン効果が生じており、且つ、不均一な外部磁界Ｂｘの影響がない場合、地磁気測定装置は、最適中心点ｃＥＯＰの示す座標、最適楕円面補正行列ＴＯＰ、及び、３次元磁気センサ６０の出力する磁気データｑｉの示す座標に基づいて、地磁気Ｂｇの向きを算出する。
このように、楕円面補正部２００、楕円面球面変換部５００、及び歪判定部９００は、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する立体ＳＤが、球面、楕円面、及び、球面とも楕円面とも異なる歪んだ形状の立体のいずれであるかを判別する、歪形状判断部４として機能し、不正確なオフセットｃＯＦＦによる不正確な地磁気Ｂｇの算出を防止することができる。
【０１７５】
また、第２実施形態に係る地磁気測定装置は、分布判定部７００を備えることにより、複数の磁気データｑ１〜ｑＮの示す座標が、センサ座標系ΣＳにおいて３次元的な広がりを有して分布しているか否かを判定することができる。これにより、複数の磁気データｑ１〜ｑＮの示す座標が２次元的または１次元的に分布している場合に、中心点算出部８００が中心点ｃＳの示す座標を算出することを防止でき、不正確な中心点ｃＳをオフセットｃＯＦＦとして採用することを防止できる。
また、複数の磁気データｑ１〜ｑＮの示す座標が２次元的または１次元的に分布している場合に、複数の磁気データｑ１〜ｑＮの示す座標の分布の形状から、楕円面の形状を特定することが困難であることが多い。従って、このような場合に、楕円面補正部２００が、楕円補正を行うことを防止する。
【０１７６】
＜Ｃ．変形例＞  
本発明は上述した実施形態に限定されるものではなく、例えば、以下の変形が可能である。また、以下に示す変形例のうちの２以上の変形例を、矛盾しない範囲で適宜組み合わせることもできる。
【０１７７】
（１）変形例１  
上述した実施形態において、歪誤差ベクトルｋ（Ｅ）に用いられる基準点ｗＫＥには、式（８９）に示した中心点ｃＳ、または式（１０７）に示した最適中心点ｃＥＯＰが採用されるが、本発明はこのような形態に限定されるものでは無く、基準点ｗＫＥとして、式（７４）に示した重心ｑＣ、または式（１０６）に示した重心ｓＣを採用してもよい。
歪評価行列Ｅの各成分は、基準点ｗＫＥの示す座標から見て、歪評価行列Ｅの最大固有値λＥ１に対応する固有ベクトルｕＥ１の方向と、球面Ｓ２からの誤差が大きな磁気データｑｉ（または、変換後磁気データｓｉ）が多く存在する領域を示す方向とが、近づくように定められる。また、歪評価行列Ｅの最大固有値λＥ１は、基準点ｗＫＥから見て、固有ベクトルｕＥ１の方向に存在する磁気データｑｉ（または、変換後磁気データｓｉ）の示す座標と球面Ｓ２との誤差の大きさを表す値となる。
従って、基準点ｗＫＥをどのような値に設定しても、複数の磁気データｑ１〜ｑＮ（または複数の変換後磁気データｓ１〜ｓＮ）の示す座標が、基準点ｗＫＥから見て広く分布している場合には、歪評価行列Ｅによって、立体ＳＤ（または立体ＳＤＥ）と球面Ｓ２との形状の相違の程度を評価することが可能となる。
【０１７８】
（２）変形例２  
上述した実施形態及び変形例において、第１の歪判定処理及び第２の歪判定処理は、いずれも、１つの基準点ｗＫＥを用いた歪誤差ベクトルｋ（Ｅ）に基づいて算出される歪評価値ｇＤ（Ｅ）により、立体ＳＤ（または立体ＳＤＥ）と球面との形状の相違の程度を評価するが、本発明はこのような形態に限定されるものではなく、２つの基準点ｗＫＥを用いて算出される異なる２つの歪誤差ベクトルｋ（Ｅ）により、２つの歪評価値ｇＤ（Ｅ）を算出して、立体ＳＤ（または立体ＳＤＥ）と球面との形状の相違の程度を評価してもよい。
例えば、第２の歪判定処理において、まず、最適中心点ｃＥＯＰを基準点ｗＫＥに採用して算出された歪評価値ｇＤ（Ｅ）に基づいて立体ＳＤＥと球面との形状の相違の程度を評価し、その後、重心ｓＣを基準点ｗＫＥに採用して算出された歪評価値ｇＤ（Ｅ）に基づいて立体ＳＤＥと球面との形状の相違の程度を評価してもよい。この場合、２度の評価の両方において、立体ＳＤＥの形状が球面の形状と看做すことができる場合に、歪判定処理の結果を肯定としてもよい。
このように、２つの基準点ｗＫＥを用いて、複数の磁気データｑ１〜ｑＮ（または、複数の変換後磁気データｓ１〜ｓＮ）の示す座標と球面Ｓ２との誤差の大きさを評価することにより、１つの基準点ｗＫＥのみを用いる場合に比べて、立体ＳＤ（または、立体ＳＤＥ）と球面との形状の相違の程度を正確に評価することができる。
【０１７９】
（３）変形例３  
上述した実施形態及び変形例において、地磁気測定装置は、第１の歪判定処理（ステップＳ５）及び第２の歪判定処理（ステップＳ９）の双方を行うが、本発明はこのような形態に限定されるものでは無く、第１の歪判定処理または第２の歪判定処理のうち、一方のみを実施してもよい。
例えば、地磁気測定装置が、第１の歪判定処理のみを行う場合、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する立体ＳＤの形状が、球面と看做すことができるか否かを判定することができるため、ソフトアイアン効果が生じているか否かを判定することができる。そして、ソフトアイアン効果が生じていない場合には、楕円面補正部２００による楕円補正を行うこと無く、３次元磁気センサ６０が出力する磁気データｑｉの示す座標と、中心点算出部８００の算出する球面Ｓの中心点ｃＳの示す座標とにより地磁気Ｂｇの向きを算出することができるため、計算負荷の低減が可能となる。
【０１８０】
また、例えば、地磁気測定装置が、第２の歪判定処理のみを行う場合、複数の変換後磁気データｓ１〜ｓＮの示す座標を近傍に有する立体ＳＤＥの形状が、球面と看做すことができるかを判定することができる。
そして、第２の歪判定処理における判定結果が肯定である場合、複数の磁気データｑ１〜ｑＮの示す座標を近傍に有する立体ＳＤの形状は、楕円面と看做すことができるため、３次元磁気センサ６０が出力する磁気データｑｉの示す座標、最適楕円面補正行列ＴＯＰ、及び、最適中心点ｃＥＯＰに基づいて、地磁気Ｂｇの向きを算出することができる。なお、楕円面には球面が含まれるため（例えば、楕円面補正行列ＴＥの有する３つの固有値が全て１の場合、等）、ソフトアイアン効果の発生の有無にかかわらず、地磁気Ｂｇの向きを算出することができる。
一方、第２の歪判定処理における判定結果が否定である場合、複数の磁気データｑ１〜ｑＮは不均一な外部磁界Ｂｘの影響を受けているため、地磁気測定装置は、オフセットｃＯＦＦ及び地磁気Ｂｇの向きの算出を防止する。
【０１８１】
（４）変形例４    
上述した実施形態及び変形例において、地磁気測定装置は、楕円面最適補正値生成部４００を備え、最適楕円面補正行列ＴＯＰ及び最適中心点ｃＥＯＰに基づいて、磁気データｑｉの示す座標を、変換後磁気データｓｉの示す座標へと変換する楕円補正を行うが、本発明はこのような形態に限定されるものでは無く、地磁気測定装置が、楕円面最適補正値生成部４００を備えずに構成されてもよい。この場合、地磁気測定装置（地磁気算出部６００または地磁気算出部６００ａ）は、楕円面初期補正値生成部３００が生成する初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０の座標を、楕円面補正行列ＴＥ及びオフセットｃＯＦＦとして採用して、楕円補正を行ってもよい。この場合、地磁気測定装置は、最適楕円面生成処理を行わないため、地磁気Ｂｇの算出に係る計算負荷を低減することができる。
なお、楕円面初期補正値生成部３００は、空間Ω上の第１評価軸ξ１、第２評価軸ξ２、及び、第３評価軸ξ３の３つの異なる評価軸を用いて、３つの異なる楕円面を生成した後に、これらの３つの楕円面の中心点の相互間の距離が第１閾値Δｃ以下であるか否かを判定し、判定結果が肯定となる場合に限り、初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０を生成する。従って、初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０により表される初期楕円面ＶＥ０は、複数の磁気データｑ１〜ｑＮの示す座標との誤差を最小化する楕円面ではないが、複数の磁気データｑ１〜ｑＮの示す座標の分布形状を正確に表す。すなわち、初期楕円面補正行列Ｔ０及び初期中心点ｃＥ０を用いた楕円補正により、正確な地磁気Ｂｇの向きを算出することができる。
【０１８２】
（５）変形例５    
上述した実施形態及び変形例において、地磁気測定装置は、第２の歪判定処理に用いられる複数の入力座標として、複数の変換後磁気データｓ１〜ｓＮの示す座標を適用したが、本発明はこのような形態に限定されるものでは無く、複数の変換後磁気データｓ１〜ｓＮの示す座標を最適中心点ｃＥＯＰを起点として表した複数のベクトル、すなわち、ベクトル（ｓ１−ｃＥＯＰ）〜ベクトル（ｓＮ−ｃＥＯＰ）を複数の入力座標に適用してもよい。この場合、歪判定処理に用いられるデータのサイズを小さくすることができ、処理に必要なメモリのサイズの節約をすると共に、処理速度を向上させるという利点を有する。
【０１８３】
（６）変形例６    
上述した実施形態及び変形例において、初期楕円面生成部３１０は、３つの楕円面（第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚ）の各々の係数行列及び中心点の座標を算出するものであったが、本発明はこのような形態に限定されるものでは無く、初期楕円面生成部３１０は、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚのうち、２つの楕円面の各々の係数行列及び中心点の座標を算出するものであってもよい。この場合、初期楕円面生成部３１０は、第１楕円面生成部３１１、第２楕円面生成部３１２、及び第３楕円面生成部３１３のうち、少なくとも２つを備えればよい。
図１３において説明したように、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚのうち２つの楕円面の形状の相違の程度（具体的には、２つの楕円面の２つの中心点間の距離）を評価することにより、複数の磁気データｑ１〜ｑＮの示す座標から楕円面の形状を特定することが困難な場合であるか否かを把握することが可能である。従って、初期楕円面生成部３１０が、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚのうち、少なくとも２つの楕円面各々の係数行列及び中心点の座標を算出すれば、不適切な初期楕円面補正行列Ｔ０が生成されることを防止することができる。また、初期楕円面生成部３１０が２つの楕円面各々の係数行列及び中心点の座標を算出する場合、３つの楕円面各々の係数行列及び中心点の座標を算出する場合に比べて、計算負荷の軽減が可能となる。
なお、初期楕円面生成部３１０が２つの楕円面の各々の係数行列及び中心点の座標を算出する場合、初期補正値生成部３３０は、当該２つの楕円面のうち、少なくとも１つの楕円面の係数行列に基づいて初期楕円面補正行列Ｔ０を算出すればよい。同様に、初期補正値生成部３３０は、当該２つの楕円面のうち、少なくとも１つの楕円面の中心点に基づいて、初期中心点ｃＥ０の座標を算出すればよい。
【０１８４】
（７）変形例７    
上述した実施形態及び変形例において、初期楕円面中心点判定部３２２は、中心点ｃｘｘ、中心点ｃｙｙ、及び中心点ｃｚｚの３つの中心点の相互間の距離の全てが、第１閾値Δｃ以下であるか否か（第２の条件を充足するか否か）を判定するものであるが、本発明はこのような判定の方法に限定されるものでは無く、初期楕円面中心点判定部３２２は、中心点ｃｘｘ、中心点ｃｙｙ、及び中心点ｃｚｚのうち２つの中心点間の距離が、第１閾値Δｃ以下であるか否かを判定するものであってもよい。
例えば、変形例６のように、初期楕円面生成部３１０が、第１楕円面Ｖｘｘ、第２楕円面Ｖｙｙ、及び第３楕円面Ｖｚｚのうち、２つの楕円面（例えは、第１楕円面Ｖｘｘ及び第２楕円面Ｖｙｙ）の各々の中心点（例えば、中心点ｃｘｘ及び中心点ｃｙｙ）の座標を算出する場合、初期楕円面中心点判定部３２２は、当該２つの中心点ｃｘｘ及びｃｙｙの間の距離が、第１閾値Δｃ以下であるか否かを判定するものであってもよい。
このような判定によっても、複数の磁気データｑ１〜ｑＮの示す座標の分布の形状から楕円面の形状を特定することが困難な場合であるか否かを把握することが可能であり、不適切な初期楕円面補正行列Ｔ０が生成されることを防止することができる。
【符号の説明】
【０１８５】
１…機器、６０…３次元磁気センサ、Ｂｇ…地磁気、Ｂｉ…内部磁界、Ｂｘ…外部磁界、Ｂｍ…着磁磁界、ｑｉ…磁気データ、ｓｉ…変換後磁気データ、２１…軟磁性材料、２…部品、ＴＥ…楕円面補正行列、ｃＯＦＦ…オフセット、２００…楕円面補正部、３００…楕円面初期補正値生成部、４００…楕円面最適補正値生成部、Ｔ０…初期楕円面補正行列、ｃＥ０…初期中心点、ＴＯＰ…最適楕円面補正行列、ｃＥＯＰ…最適中心点、３１０…初期楕円面生成部、３２０…初期楕円面判定部、３３０…初期補正値生成部、ｃ…変数ベクトル、Ｔ…変数行列、８００…中心点算出部、９００…歪判定部、ｃＳ…中心点。
【図１】
0005866893.tif 000039
【図２】
0005866893.tif 000040
【図３】
0005866893.tif 000041
【図４】
0005866893.tif 000042
【図５】
0005866893.tif 000043
【図６】
0005866893.tif 000044
【図７】
0005866893.tif 000045
【図８】
0005866893.tif 000046
【図９】
0005866893.tif 000047
【図１０】
0005866893.tif 000048
【図１１】
0005866893.tif 000049
【図１２】
0005866893.tif 000050
【図１３】
0005866893.tif 000051
【図１４】
0005866893.tif 000052
【図１５】
0005866893.tif 000053
【図１６】
0005866893.tif 000054
【図１７】
0005866893.tif 000055
【図１８】
0005866893.tif 000056
【図１９】
0005866893.tif 000057
【図２０】
0005866893.tif 000058
【図２１】
0005866893.tif 000059
【図２２】
0005866893.tif 000060
【図２３】
0005866893.tif 000061
【図２４】
0005866893.tif 000062
【図２５】
0005866893.tif 000063