JP5266846B2 - 測位方法、プログラム及び測位装置 - Google Patents

Description

本発明は、測位装置が、各カルマンフィルタ処理それぞれをモデルとみなし、前記カルマンフィルタ処理それぞれの処理結果を所定のモデル確率で重みづけして合成するＩＭＭ演算を行って測位する測位方法等に関する。

人工衛星を利用した測位システムとしては、ＧＰＳ（Global Positioning System）が広く知られており、携帯型電話機やカーナビゲーション装置等に内蔵された測位装置に利用されている。ＧＰＳでは、自機の位置を示す３次元の座標値と、時計誤差との４つのパラメータの値を、複数のＧＰＳ衛星の位置や各ＧＰＳ衛星から自機までの擬似距離等の情報に基づいて求める測位演算を行うことで、自機の現在位置を測位する。

測位演算としては、最小二乗法を用いた測位演算や、カルマンフィルタを用いた測位演算が知られている。そして、カルマンフィルタを用いた測位演算を応用した技術として、測位装置の移動状態に対応する複数のカルマンフィルタモデルを用意して各モデルについて測位演算を行い、その測位演算の結果を合成して測位位置を決定するインタラクティブミキシングモデル（Interactive Mixing Model）演算（以下、「ＩＭＭ演算」と称す。）が知られている（例えば、非特許文献１）。
Yaakov Bar-Shalom, X.Rong Li, Thiagalingam Kirubarajan著、「Estimation with Applications to Tracking and Navigation」、（米国）、JOHN WILEY ＆ SONS INC、２００１年、p.９２−９３、４５３−４５７

ＩＭＭ演算による測位では、カルマンフィルタの入力値となる外部観測量の尤もらしさ（以下、「尤度」と称す。）を、尤度関数と呼ばれる正規分布（ガウス分布）に従った関数から算出し、この算出した尤度を用いて、各カルマンフィルタモデルの確率（以下、「モデル確率」と称す。）を算出する。そして、算出したモデル確率を用いて各カルマンフィルタモデルの処理結果を合成することで、測位装置の現在位置を求める。

しかし、尤度関数はＥＸＰ（Exponential）関数で表されるため、ＥＸＰ関数の指数部分が非常に小さな値（例えば「−２０以下」）になると、尤度は限りなく「０」に近くなり、計算機で計算しようとしても演算可能な数値範囲の関係から結局「０」になってしまう。モデル確率は尤度を用いて算出するため、尤度が「０」になってしまうと、モデル確率を正しく算出することができず、その結果、測位位置を求めることができないという大きな問題があった。

本発明は、上述した課題に鑑みて為されたものである。

以上の課題を解決するための第１の発明は、測位装置が、各カルマンフィルタ処理それぞれをモデルとみなし、前記カルマンフィルタ処理それぞれの処理結果を所定のモデル確率で重みづけして合成するＩＭＭ演算を行って測位する測位方法であって、前記カルマンフィルタ処理それぞれの尤度指標値を算出することと、前記算出された尤度指標値の前記カルマンフィルタ処理間の相対値を算出することと、前記算出された相対値を用いて前記カルマンフィルタ処理それぞれの前記モデル確率を算出することと、前記算出されたモデル確率を用いて前記カルマンフィルタ処理それぞれの処理結果を合成することで測位位置を決定することと、を含む測位方法である。

また、他の発明として、各カルマンフィルタ処理それぞれをモデルとみなし、前記カルマンフィルタ処理それぞれの処理結果を所定のモデル確率で重みづけして合成するＩＭＭ演算を行って測位する測位装置であって、前記カルマンフィルタ処理それぞれの尤度指標値を算出する指標値算出部と、前記算出された尤度指標値の前記カルマンフィルタ処理間の相対値を算出する相対値算出部と、前記算出された相対値を用いて前記カルマンフィルタ処理それぞれの前記モデル確率を算出するモデル確率算出部と、前記算出されたモデル確率を用いて前記カルマンフィルタ処理それぞれの処理結果を合成することで測位位置を決定する測位位置決定部と、を備えた測位装置を構成してもよい。

この第１の発明等によれば、カルマンフィルタ処理それぞれの尤度指標値が算出され、この尤度指標値のカルマンフィルタ処理間の相対値が算出される。そして、算出された相対値を用いてカルマンフィルタ処理それぞれのモデル確率が算出され、このモデル確率を用いてカルマンフィルタ処理それぞれの処理結果が合成されることで測位位置が決定される。

ここで、尤度指標値とは、カルマンフィルタの入力値となる外部観測量の尤度を算出するための尤度関数の指数部分のことを指す。ＩＭＭ演算におけるモデル確率は、尤度指標値を用いて算出した外部観測量の尤度を基に算出することができるが、尤度指標値が非常に小さな値である場合には、計算機で計算しようとしても演算可能な数値範囲の関係で尤度が「０」になってしまい、モデル確率を正しく求めることができない。

本願発明者は、この問題を解決するため、尤度指標値をそのまま用いるのではなく、複数のカルマンフィルタ処理間の尤度指標値の相対値を用いることで、モデル確率を正しく求めることができることを発見した。これは、個々の尤度指標値が非常に小さな値であったとしても、尤度指標値の相対値は概ね「０」の近傍の値となり、指数関数に当該相対値を与えて計算すると「０」とはならない点に着目したものである。かかる原理により、モデル確率を正しく算出することが可能となり、その結果、正確な測位位置の算出を実現し得る。

また、第２の発明として、第１の発明の測位方法であって、前記モデル確率を算出することは、前記相対値が小さいほど値が小さくなる所定の指数関数を用いて前記モデル確率を算出することである測位方法を構成してもよい。

この第２の発明によれば、カルマンフィルタ処理間の尤度指標値の相対値が小さいほど値が小さくなる指数関数（例えば、尤度指標値の相対値を変数とするＥＸＰ関数）を用いてモデル確率が算出される。

また、第３の発明として、第２の発明の測位方法であって、前記モデル確率を算出することは、前記相対値を用いて算出した前記所定の指数関数の値が極限値相当値を示す条件である所定の極限値条件を満たすか否かを判定することと、前記所定の極限値条件を満たす場合に、前記指数関数の値を極限値とみなして前記モデル確率を算出することと、を含む測位方法を構成してもよい。

この第３の発明によれば、カルマンフィルタ処理間の尤度指標値の相対値を用いて算出した指数関数の値が極限値条件を満たす場合は、当該指数関数の値を極限値とみなしてモデル確率が算出される。これにより、指数関数の値が０に収束したり、＋∞に発散する場合であっても、計算機によってモデル確率を求めることが可能となる。

また、第４の発明として、第３の発明の測位方法であって、前記所定の極限値条件には、最大値相当値を示す条件と最小値相当値を示す条件とが含まれ、前記モデル確率を算出することは、前記所定の極限値条件を満たす場合に、前記最大値相当値を示す条件を満たすか、前記最小値相当値を示す条件を満たすかに応じて、前記モデル確率を１又は０とすることを含む測位方法を構成してもよい。

この第４の発明によれば、カルマンフィルタ処理間の尤度指標値の相対値を用いて算出した指数関数の値が最大値相当値を示す条件を満たすか、最小値相当値を満たすかに応じて、モデル確率が１又は０とされる。

また、第５の発明として、第１〜第４の何れかの発明の測位方法であって、前記複数のカルマンフィルタ処理には、前記測位装置の移動状態として、停止状態にあることを前提としたカルマンフィルタ処理と、等速移動状態にあることを前提としたカルマンフィルタ処理と、定加速度状態にあることを前提としたカルマンフィルタ処理との３種類のカルマンフィルタ処理が少なくとも含まれる測位方法を構成してもよい。

この第５の発明によれば、測位装置の移動状態として、停止状態と、等速移動状態と、定加速度状態とに対応する３種類のカルマンフィルタ処理が行われ、これらのカルマンフィルタ処理の処理結果が合成される。

また、第６の発明として、第１〜第５の何れかの発明の測位方法を前記測位装置に内蔵されたコンピュータに実行させるためのプログラムを構成してもよい。

１．原理
図１は、本実施形態におけるインタラクティブミキシングモデル（Interactive Mixing Model：以下、「ＩＭＭ」と称す。）の概要を説明するための図である。ここでは、説明を分かりやすくするため、複数の処理回路ブロックにおける演算によってＩＭＭ演算が実現されるものとして説明するが、単独のプロセッサによってＩＭＭ演算を実現することも可能であり、図２及び図３の処理フローの説明においては単独のプロセッサにより実行するものとして説明する。

ＩＭＭは、複数のカルマンフィルタモデル（Kalman Filter Model：以下、「ＫＦモデル」と称す。）それぞれの処理結果をモデル確率「μ」と呼ばれる重みで合成することで、最終的な出力値を決定するモデルである。図１では、第１ＫＦモデル〜第ｎＫＦモデルの「ｎ個」のＫＦモデルを含むＩＭＭを図示している。

また、説明の簡明化のため、各ＫＦモデルに対して、単一の外部観測量（本実施形態では「メジャメント実測値」が相当する。）「Ｚ」が与えられるものとして説明するが、複数のメジャメント実測値が与えられる場合についても同様である。メジャメント実測値は、ＩＭＭ演算を適用するシステムに応じて適宜設定することができる。

ＩＭＭ演算では、ミキサが、各ＫＦモデル間の遷移確率を示すモデル間遷移確率「ｐ」と、１つ前のステップにおいてモデル確率更新部により更新されたモデル確率「μ：μ₁，μ₂・・・，μ_n」と、各ＫＦモデルそれぞれにおいて算出された状態ベクトルの計算値「Ｘ：Ｘ₁，Ｘ₂，・・・，Ｘ_n」及び誤差共分散の計算値「Ｐ：Ｐ₁，Ｐ₂，・・・，Ｐ_n」とを用いて、各ＫＦモデルそれぞれの状態ベクトル及び誤差共分散の今回の計算における初期値「Ｘ₀：Ｘ₀₁，Ｘ₀₂，・・・，Ｘ_0n」及び「Ｐ₀：Ｐ₀₁，Ｐ₀₂，・・・，Ｐ_0n」を算出し、計算結果を、対応するＫＦモデルにそれぞれ出力する。

また、ミキサは、モデル間遷移確率「ｐ」と、１つ前のステップにおいてモデル確率更新部により更新されたモデル確率「μ」とを用いて、規格化定数「ＣＮ：Ｃ₁，Ｃ₂，・・・，Ｃ_n」と呼ばれる定数を算出して、モデル確率更新部に出力する。

各ＫＦモデルは、ミキサから入力した初期値「Ｘ₀」及び「Ｐ₀」と、メジャメント実測値「Ｚ：Ｚ₁，Ｚ₂，・・・，Ｚ_n」とを用いて、公知のカルマンフィルタと同様、予測演算と補正演算とを行うことで、状態ベクトル及び誤差共分散の今回の計算値を算出して、モデル確率更新部に出力する。また、各ＫＦモデルは、尤度指標値「Ｖｘ：Ｖｘ₁，Ｖｘ₂，・・・，Ｖｘ_n」と重み係数「β：β₁，β₂，・・・，β_n」と呼ばれる値を算出して、モデル確率更新部に出力する。

尤度指標値「Ｖｘ」は、メジャメント実測値「Ｚ」の尤もらしさ（以下、「尤度」と称す。）を算出するための関数である尤度関数の指数部分として表される値である。また、重み係数「β」は、尤度関数の係数部分として表される値である。

モデル確率更新部は、ミキサから入力した規格化定数「ＣＮ」と、ＫＦモデルから入力した尤度指標値「Ｖｘ」及び重み係数「β」とを用いて、モデル確率「μ」を新たに算出・更新する。そして、ＫＦモデルから入力した状態ベクトル及び誤差共分散それぞれを、算出したモデル確率「μ」で重み付けして合成することで、状態ベクトルの出力値「Ｘ」及び誤差共分散の出力値「Ｐ」を算出する。

次に、実際の数式を参照しながら、単一のプロセッサでＩＭＭ演算を実現する場合の処理の流れについて説明する。以下の数式において、各変数の添え字の「ｊ」は、ｊ番目のＫＦモデルに対応する変数であることを示しているものとする。

図２は、プロセッサの一種であるＣＰＵ（Central Processing Unit）が実行するＩＭＭ演算処理の流れを示すフローチャートである。
先ず、ＣＰＵは、モデル間遷移確率「ｐ」を取得する（ステップＳ１）。モデル間遷移確率「ｐ」は、次式（１）に示すようなｎ×ｎの行列で与えられる。

行列の成分「ｐ_nm」は、ｎ番目のＫＦモデルからｍ番目のＫＦモデルへの遷移の確率を示している。例えば、行列の成分「ｐ₁₁」は、第１ＫＦモデルから第１ＫＦモデル（自身）への遷移の確率を示しており、行列の成分「ｐ₁₂」は、第１ＫＦモデルから第２ＫＦモデルへの遷移の確率を示している。また、各行「ｉ」の成分を合算した値は「１」である。すなわち、「ｐ_i1＋ｐ_i2＋・・・＋ｐ_in＝１」である。

次いで、ＣＰＵは、次式（２）に従って規格化定数「ＣＮ」を算出する（ステップＳ３）。

式（２）からわかるように、規格化定数「ＣＮ」は、モデル間遷移確率「ｐ」の転置行列と、モデル確率「μ」との積で表される。

次いで、ＣＰＵは、次式（３）に従って、合成確率「γ」を算出する（ステップＳ５）。

式（３）からわかるように、合成確率「γ」は、モデル間遷移確率「ｐ」の各成分と、モデル確率「μ」の各成分と、規格化定数「ＣＮ」の各成分とを用いて算出される。

その後、ＣＰＵは、次式（４）及び（５）に従って、ＫＦ処理における状態ベクトルの初期値「Ｘ₀」及び誤差共分散の初期値「Ｐ₀」を算出する（ステップＳ７）。

次いで、ＣＰＵは、各ＫＦモデルそれぞれについて、ループＡの処理を実行する（ステップＳ９〜Ｓ２３）。ループＡでは、ＣＰＵは、次式（６）及び（７）に従って予測演算を行って、状態ベクトル及び誤差共分散を予測する（ステップＳ１１）。

但し、添え字の「−」は予測値であることを示している。また、「φ」及び「Ｑ」は、それぞれ状態遷移行列及びプロセスノイズと呼ばれる行列であり、各ＫＦモデル毎に設定される行列である。

次いで、ＣＰＵは、外部（例えば外部観測量を検出するセンサ）からメジャメント実測値「Ｚ」を取得する（ステップＳ１３）。また、ＣＰＵは、ステップＳ１１で予測した状態ベクトルの予測値「Ｘ⁻」を用いて、外部観測量の予測値（以下、「メジャメント予測値」と称す。）を算出する（ステップＳ１５）。メジャメント予測値は、所定の観測行列「Ｈ」と状態ベクトルの予測値「Ｘ⁻」との積「ＨＸ⁻」で表される。

その後、ＣＰＵは、補正演算を行って、ステップＳ１１で予測した状態ベクトル及び誤差共分散を補正する（ステップＳ１７）。具体的には、先ず次式（８）に従って、カルマンゲイン「Ｋ」と呼ばれるゲインを算出する。

但し、Ｅ［・］は平均値（期待値）を示している。また、「Ｒ」は測定誤差行列と呼ばれるメジャメント実測値「Ｚ」に含まれる誤差を示す行列である。

そして、ＣＰＵは、メジャメント実測値「Ｚ」とメジャメント予測値「ＨＸ⁻」とを用いて、次式（９）及び（１０）に従って、状態ベクトル及び誤差共分散を補正する。

但し、添え字の「＋」は補正値であることを示している。また、「Ｉ」は単位行列である。

次いで、ＣＰＵは、次式（１１）に従って、尤度指標値「Ｖｘ」を算出する。

また、ＣＰＵは、次式（１２）に従って、重み係数「β」を算出する。

但し、「Ｎ」は取得したメジャメント実測値の個数を示している。

全てのＫＦモデルについてステップＳ１１〜Ｓ２１の処理を行った後、ＣＰＵは、ループＡを終了する。その後、ＣＰＵは、各ＫＦモデル同士の相対的な組合せそれぞれにおける尤度指標値「Ｖｘ」の相対値を算出する（ステップＳ２５）。例えば、第１ＫＦモデル〜第３ＫＦモデルの３つのＫＦモデルが存在する場合は、「Ｖｘ₁−Ｖｘ₂」、「Ｖｘ₂−Ｖｘ₁」、「Ｖｘ₁−Ｖｘ₃」、「Ｖｘ₃−Ｖｘ₁」、「Ｖｘ₂−Ｖｘ₃」及び「Ｖｘ₃−Ｖｘ₂」の６つの相対的な組合せに対して、それぞれの値を、尤度指標値の相対値として算出する。

その後、ＣＰＵは、ステップＳ２５で算出した相対値それぞれが、所定の最大値相当値（例えば「＋２０」）以上となったか、また、所定の最小値相当値（例えば「−２０」）以下となったか否かを判定する（ステップＳ２７）。

そして、ＣＰＵは、次式（１３−２）で与えられるモデル確率「μ」の式の分母において、相対値が最大値相当値以上となったと判定した項は「＋∞」、相対値が最小値相当値以下となったと判定した項は「０」として、各ＫＦモデルそれぞれについてモデル確率「μ」を算出する（ステップＳ２９）。

より具体的には、相対値が最大値相当値以上となった項が１つでも存在する場合は、モデル確率「μ」を「０」とし、全ての項の相対値が最小値相当値以下となった場合は、モデル確率「μ」を「１」とする。

尚、式（１３−１）における「λ」は、次式（１４）で与えられる尤度であり、「Ｃ」は、次式（１５）で与えられる規格化定数である。

次いで、ＣＰＵは、ステップＳ２９で算出したモデル確率「μ」を用いて、次式（１６）及び（１７）に従って、状態ベクトル及び誤差共分散の出力値「Ｘ」及び「Ｐ」を決定する（ステップＳ３１）。そして、ＣＰＵは、ＩＭＭ演算処理を終了する。

２．実験結果
次に、ＫＦモデルが、第１ＫＦモデル及び第２ＫＦモデルの２つである場合に、各ＫＦモデルそれぞれについてモデル確率「μ」を算出した場合の実験結果について説明する。ここでは、第１ＫＦモデルの尤度指標値「Ｖｘ₁」を「−２０」、第１ＫＦモデルの尤度指標値「Ｖｘ₂」を「−２１」として実験を行った。ＫＦモデルが２つである場合は、次式（１８）及び（１９）に従ってモデル確率「μ₁」及び「μ₂」を算出する。

図４は、従来の手法に則って、式（１８−１）及び（１９―１）によりモデル確率「μ₁」及び「μ₂」を算出した場合の各種パラメータの値を示した図である。この結果を見ると、「λ₁」及び「λ₂」が共に「０」であるために、「λ₁Ｃ₁」及び「λ₂Ｃ₂」が「０」となり、「Ｃ＝λ₁Ｃ₁＋λ₂Ｃ₂」も「０」となっている。これにより、式（１８−１）及び（１９−１）の分母がそれぞれ「０」となるため、モデル確率「μ₁」及び「μ₂」は不定となっている。

この様子を定性的に示したグラフが、図６である。同図において、横軸は「Ｖｘ₁」、縦軸は「ｙ＝ｅｘｐ（Ｖｘ₁）」をそれぞれ示している。上記の実験では、「Ｖｘ₁」の値が「−２０」であるが、このグラフを見ると、「ｙ」の値はほぼ「０」となることがわかる。また、「Ｖｘ₂」の値は「−２１」であるが、この場合も同様である。

すなわち、従来の手法では、各ＫＦモデルそれぞれについて算出された尤度指標値「Ｖｘ」が非常に小さな負の領域（図６のハッチング領域）に集中すると、ＥＸＰ関数の値が演算可能な数値範囲の関係で「０」となり、式（１４）により算出される尤度「λ」は「０」になってしまうという問題がある。尤度が「０」になると、式（１３−１）によってモデル確率「μ」を算出することができない。

図５は、本実施形態の手法に則って、式（１８−２）及び（１９−２）によりモデル確率「μ₁」及び「μ₂」を算出した場合の各種パラメータの値を示した図である。この結果を見ると、尤度指標値「Ｖｘ₁」及び「Ｖｘ₂」ではなく、その相対値「Ｖｘ₂−Ｖｘ₁」及び「Ｖｘ₁−Ｖｘ₂」を用いて計算を行ったことで、モデル確率「μ₁」及び「μ₂」が正しく求められていることがわかる。

この様子を定性的に示したグラフが、図７である。同図において、横軸は「Ｖｘ₁−Ｖｘ₂」、縦軸は「ｙ＝ｅｘｐ（Ｖｘ₁−Ｖｘ₂）」をそれぞれ示している。上記の実験では、「Ｖｘ₁」の値が「−２０」であり、「Ｖｘ₂」の値が「−２１」であるため、「Ｖｘ₁−Ｖｘ₂」の値は「１」となるが、この場合は「ｙ」の値を正しく求めることができる。また、「Ｖｘ₂−Ｖｘ₁」の値は「−１」となるが、この場合も同様に「ｙ」の値を正しく求めることができる。

本願発明者は、各ＫＦモデルそれぞれについて算出された尤度指標値「Ｖｘ」が非常に小さな値であったとしても、尤度指標値の相対値は概ね「０」の近傍の領域（図７のハッチング領域）に集中する点に着目した。そして、式（１４）に従って尤度「λ」を算出することなく、式（１３−１）を展開・変形することで得られる式（１３−２）を用いて、尤度指標値の相対値からモデル確率「μ」を直接算出することで、尤度指標値「Ｖｘ」の大きさに左右されずにモデル確率「μ」を求めることができることを発見した。

３．実施例
次に、測位装置を備えた電子機器の一種である携帯型電話機に上述した原理を適用した場合の実施例について説明する。但し、上述した原理の適用可能な実施例がこれに限定されるわけではない。

３−１．機能構成
図８は、本実施例における携帯型電話機１の機能構成を示すブロック図である。携帯型電話機１は、ＧＰＳアンテナ１０と、ＧＰＳ受信部２０と、ＴＣＸＯ（Temperature Compensated Crystal Oscillator）４０と、ホストＣＰＵ（Central Processing Unit）５０と、操作部６０と、表示部７０と、携帯電話用アンテナ８０と、携帯電話用無線通信回路部９０と、ＲＯＭ（Read Only Memory）１００と、ＲＡＭ（Random Access Memory）１１０とを備えて構成される。

ＧＰＳアンテナ１０は、ＧＰＳ衛星から発信されているＧＰＳ衛星信号を含むＲＦ（Radio Frequency）信号を受信するアンテナであり、受信した信号をＧＰＳ受信部２０に出力する。尚、ＧＰＳ衛星信号は、衛星毎に異なる拡散符号の一種であるＰＲＮ（Pseudo Random Noise）コードで直接スペクトラム拡散方式により変調された１．５７５４２［ＧＨｚ］の通信信号である。ＰＲＮコードは、コード長１０２３チップを１ＰＮフレームとする繰返し周期１ｍｓの擬似ランダム雑音符号である。

ＧＰＳ受信部２０は、ＧＰＳアンテナ１０から出力された信号に基づいて携帯型電話機１の現在位置を測位する測位回路であり、いわゆるＧＰＳ受信機に相当する機能ブロックである。ＧＰＳ受信部２０は、ＲＦ（Radio Frequency）受信回路部２１と、ベースバンド処理回路部３０とを備えて構成される。尚、ＲＦ受信回路部２１と、ベースバンド処理回路部３０とは、それぞれ別のＬＳＩ（Large Scale Integration）として製造することも、１チップとして製造することも可能である。

ＲＦ受信回路部２１は、ＲＦ信号の処理回路ブロックであり、ＴＣＸＯ４０により生成された発振信号を分周或いは逓倍することで、ＲＦ信号乗算用の発振信号を生成する。そして、生成した発振信号を、ＧＰＳアンテナ１０から出力されたＲＦ信号に乗算することで、ＲＦ信号を中間周波数の信号（以下、「ＩＦ（Intermediate Frequency）信号」と称す。）にダウンコンバートし、ＩＦ信号を増幅等した後、Ａ／Ｄ変換器でデジタル信号に変換して、ベースバンド処理回路部３０に出力する。

ベースバンド処理回路部３０は、ＲＦ受信回路部２１から出力されたＩＦ信号に対して相関処理等を行ってＧＰＳ衛星信号を捕捉・抽出し、データを復号して航法メッセージや時刻情報等を取り出して測位演算を行う回路部である。ベースバンド処理回路部３０は、演算制御部３１と、ＲＯＭ３５と、ＲＡＭ３７とを備えて構成される。また、演算制御部３１は、メジャメント取得演算部３２と、測位演算部３３とを備えて構成される。

尚、メジャメント取得演算部３２と、測位演算部３３とは、それぞれ別のＬＳＩとして製造することも、１チップとして製造することも可能である。また、本実施形態においては現在位置の測位演算そのものは測位演算部３３で実行することとして説明するが、測位演算部３３で実行する処理の一部又は全部をホストＣＰＵ５０で実行することとしてもよいのは勿論である。

メジャメント取得演算部３２は、ＲＦ受信回路部２１から出力された受信信号（ＩＦ信号）から、ＧＰＳ衛星信号の捕捉・追尾を行う回路部であり、相関演算部３２１を備えて構成されている。メジャメント取得演算部３２は、捕捉・追尾したＧＰＳ衛星信号の受信周波数やコード位相等の情報を取得し、メジャメント実測値として測位演算部３３に出力する。

相関演算部３２１は、受信信号に含まれるＰＲＮコードとレプリカコードとの相関を、例えばＦＦＴ演算を用いて算出するコヒーレント処理を行い、このコヒーレント処理の結果である相関値を所定秒数分（例えば「１秒分」）積算して積算相関値を算出するインコヒーレント処理を行うことで、ＧＰＳ衛星信号を捕捉する。レプリカコードとは、擬似的に発生させた捕捉しようとするＧＰＳ衛星信号に含まれるＰＲＮコードを模擬した信号である。

捕捉しようとするＧＰＳ衛星信号が間違いなければ、そのＧＰＳ衛星信号に含まれるＰＲＮコードとレプリカコードとは一致し（捕捉成功）、間違っていれば一致しない（捕捉失敗）。そのため、算出された積算相関値のピークを判定することによってＧＰＳ衛星信号の捕捉が成功したか否かを判定でき、レプリカコードを次々に変更して、同じ受信信号との相関演算を行うことで、ＧＰＳ衛星信号を捕捉することが可能となる。

また、相関演算部３２１は、上述したコヒーレント処理を、レプリカコードの発生信号の周波数、及び、ＰＲＮコードとレプリカコードとを相関演算する際の位相を変更しつつ行っている。レプリカコードの発生信号の周波数と受信信号の周波数とが一致し、且つ、ＰＲＮコードとレプリカコードとの位相が一致した場合に、積算相関値が最大となる。

より具体的には、捕捉対象のＧＰＳ衛星信号に応じた所定の周波数及びコード位相の範囲をサーチ範囲として設定する。そして、このサーチ範囲内で、ＰＲＮコードの開始位置（コード位相）を検出するための位相方向の相関演算と、周波数を検出するための周波数方向の相関演算とを行う。サーチ範囲は、周波数についてはＧＰＳ衛星信号の搬送波周波数である１．５７５４２［ＧＨｚ］を中心とする所定の周波数掃引範囲、コード位相についてはＰＲＮコードのチップ長である１０２３チップのコード位相範囲内に定められる。

測位演算部３３は、メジャメント取得演算部３２から入力したメジャメント実測値を用いてＩＭＭ演算を行って携帯型電話機１の現在位置を測位する測位演算を行うプロセッサである。

ＴＣＸＯ４０は、所定の発振周波数で発振信号を生成する温度補償型水晶発振器であり、生成した発振信号をＲＦ受信回路部２１及びベースバンド処理回路部３０に出力する。

ホストＣＰＵ５０は、ＲＯＭ１００に記憶されているシステムプログラム等の各種プログラムに従って携帯型電話機１の各部を統括的に制御するプロセッサである。ホストＣＰＵ５０は、測位演算部３３から入力した出力位置をプロットしたナビゲーション画面を、表示部７０に表示させる。

操作部６０は、例えばタッチパネルやボタンスイッチ等により構成される入力装置であり、押下されたキーやボタンの信号をホストＣＰＵ５０に出力する。この操作部６０の操作により、通話要求やメールの送受信要求等の各種指示入力がなされる。

表示部７０は、ＬＣＤ（Liquid Crystal Display）等により構成され、ホストＣＰＵ５０から入力される表示信号に基づいた各種表示を行う表示装置である。表示部７０には、ナビゲーション画面や時刻情報等が表示される。

携帯電話用アンテナ８０は、携帯型電話機１の通信サービス事業者が設置した無線基地局との間で携帯電話用無線信号の送受信を行うアンテナである。

携帯電話用無線通信回路部９０は、ＲＦ変換回路、ベースバンド処理回路等によって構成される携帯電話の通信回路部であり、携帯電話用無線信号の変調・復調等を行うことで、通話やメールの送受信等を実現する。

ＲＯＭ１００は、ホストＣＰＵ５０が携帯型電話機１を制御するためのシステムプログラムや、ナビゲーション機能を実現するための各種プログラムやデータ等を記憶している。

ＲＡＭ１１０は、ホストＣＰＵ５０により実行されるシステムプログラム、各種処理プログラム、各種処理の処理中データ、処理結果などを一時的に記憶するワークエリアを形成している。

３−２．データ構成
図９は、ベースバンド処理回路部３０のＲＯＭ３５に格納されたデータの一例を示す図である。ＲＯＭ３５には、測位演算部３３により読み出され、ベースバンド処理（図１３参照）として実行されるベースバンド処理プログラム３５１が記憶されている。また、ベースバンド処理プログラム３５１には、ＩＭＭ演算処理（図２及び図３参照）として実行されるＩＭＭ演算プログラム３５１１がサブルーチンとして含まれている。

ベースバンド処理とは、測位演算部３３が、ＩＭＭ演算処理を行って出力位置を決定し、当該出力位置をホストＣＰＵ５０に出力する処理である。ベースバンド処理については、フローチャートを用いて詳細に後述する。

ＩＭＭ演算処理とは、測位演算部３３が、複数のＫＦモデルに基づくＫＦ処理の処理結果を、算出・更新したモデル確率で合成することで、出力位置を決定する処理である。本実施例では、図１１に示すように、携帯型電話機１の移動状態を示すＫＦモデルとして、停止状態にあることを示す定点ＫＦモデルと、等速移動状態にあることを示す等速度ＫＦモデルと、定加速度状態にあることを示す定加速度ＫＦモデルとの３種類のＫＦモデルの処理結果に基づくＩＭＭ演算を行って現在位置を測位する。

図１０は、ベースバンド処理回路部３０のＲＡＭ３７に格納されるデータの一例を示す図である。ＲＡＭ３７には、ＩＭＭパラメータデータ３７１と、捕捉衛星別メジャメントデータ３７３と、出力位置データ３７５とが記憶される。

ＩＭＭパラメータデータ３７１は、ＩＭＭ演算処理において用いられるＩＭＭ演算の各種パラメータ（以下、「ＩＭＭパラメータ」と称す。）の値が記憶されたデータであり、ＩＭＭ演算処理において測位演算部３３により更新される。

図１２は、捕捉衛星別メジャメントデータ３７３のデータ構成例を示す図である。捕捉衛星別メジャメントデータ３７３は、ＩＭＭ演算処理に使用するメジャメント情報についてのデータであり、捕捉衛星３７３１と、メジャメント実測値３７３３と、メジャメント予測値３７３５とが対応付けて記憶される。捕捉衛星３７３１には、当該捕捉衛星の番号が記憶され、メジャメント実測値３７３３及びメジャメント予測値３７３５には、当該捕捉衛星から受信したＧＰＳ衛星信号の受信周波数やコード位相の実測値及び予測値がそれぞれ記憶される。

例えば、捕捉衛星「Ｓ１」についてのメジャメント実測値は、受信周波数が「ＳＦｒｅｑ１」、コード位相が「ＳＣＰ１」であり、メジャメント予測値は、受信周波数が「ＥＦｒｅｑ１」、コード位相が「ＥＣＰ１」である。

ＩＭＭ演算処理では、測位演算部３３は、メジャメント取得演算部３２により取得演算されたＧＰＳ衛星信号の受信周波数及びコード位相の情報を、メジャメント実測値として取得する。また、測位演算部３３は、ＫＦ処理において予測演算により予測した状態ベクトルの予測値と観測行列とを用いてＧＰＳ衛星信号の受信周波数及びコード位相の情報を予測して、メジャメント予測値とする。

出力位置データ３７５は、ホストＣＰＵ５０に出力する位置として決定された出力位置のデータであり、ＩＭＭ演算処理において測位演算部３３により更新される。

３−３．処理の流れ
図１３は、測位演算部３３によりＲＯＭ３５に記憶されているベースバンド処理プログラム３５１が読み出されて実行されることで、携帯型電話機１において実行されるベースバンド処理の流れを示すフローチャートである。ベースバンド処理は、ＲＦ受信回路部２１によるＧＰＳ衛星信号の受信と併せて、測位演算部３３が、操作部６０に測位開始指示の操作がなされたことを検出した場合に実行を開始する処理であり、各種アプリケーションの実行といった各種の処理と並行して行われる処理である。

尚、携帯型電話機１の電源のＯＮ／ＯＦＦとＧＰＳの起動／停止とを連動させ、携帯型電話機１の電源投入操作を検出した場合に処理の実行を開始させることにしてもよい。原則として、測位演算は「１秒」毎に行われるものとする。また、特に説明しないが、以下のベースバンド処理の実行中は、ＧＰＳアンテナ１０によるＲＦ信号の受信や、ＲＦ受信回路部２１によるＩＦ信号へのダウンコンバート、メジャメント取得演算部３２によるＧＰＳ衛星信号の受信周波数やコード位相の情報の取得・算出等が随時行われている状態にあるものとする。

先ず、測位演算部３３は、初期設定を行う（ステップＡ１）。具体的には、ＩＭＭパラメータを初期化し、ＲＡＭ３７のＩＭＭパラメータデータ３７１に記憶させる。次いで、測位演算部３３は、ＲＯＭ３５に記憶されているＩＭＭ演算プログラム３５１１を読み出して実行することで、ＩＭＭ演算処理を行う（ステップＡ３）。

ステップＡ３のＩＭＭ演算処理では、測位演算部３３は、図２及び図３のフローチャートに従って処理を実行する。この場合において、測位演算部３３は、携帯型電話機１の３次元の位置「（ｘ，ｙ，ｚ）」、３次元の速度「（ｕ，ｖ，ｗ）」、３次元の加速度「（ａ_x，ａ_y，ａ_z）」、クロックバイアス「ｂ」及びクロックドリフト「ｄ」を成分とする次式（２０）で表される１１次元の状態ベクトル「Ｘ」を用いてＩＭＭ演算を行う。

また、測位演算部３３は、定点ＫＦモデル用の状態遷移行列として次式（２１）で与えられる「φ₁」、等速度ＫＦモデル用の状態遷移行列として次式（２２）で与えられる「φ₂」、定加速度ＫＦモデル用の状態遷移行列として次式（２３）で与えられる「φ₃」を用いてＫＦ処理を行う。

より具体的には、測位演算部３３は、図２のＩＭＭ演算処理のステップＳ１１における予測演算において、状態遷移行列「φ₁」〜「φ₃」を用いて、式（６）及び式（７）に従って状態ベクトル及び誤差共分散の予測値「Ｘ⁻：Ｘ⁻ ₁，Ｘ⁻ ₂，Ｘ⁻ ₃」及び「Ｐ⁻：Ｐ⁻ ₁，Ｐ⁻ ₂，Ｐ⁻ ₃」を算出する。尚、状態遷移行列は、１１×１１の行列であり、その行と列の並びは、式（２０）の状態ベクトル「Ｘ」の成分の並びに対応している。

式（２１）の状態遷移行列「φ₁」では、３次元の位置に対応する３×３の行列部分の対角成分は「１」になっているが、３次元の速度及び加速度に対応する６×６の行列部分の成分は全て「０」になっており、携帯型電話機１が停止状態にあることを表していることがわかる。

式（２２）の状態遷移行列「φ₂」では、３次元の位置及び速度に対応する６×６の行列部分の対角成分は「１」になっているが、３次元の加速度に対応する３×３の行列部分の成分は全て「０」になっており、携帯型電話機１が等速移動状態にあることを表していることがわかる。

式（２３）の状態遷移行列「φ₃」では、３次元の位置、速度及び加速度に対応する９×９の行列部分の対角成分が「１」となっており、携帯型電話機１が定加速度状態にあることを表していることがわかる。

また、測位演算部３３は、定点ＫＦモデル、等速度ＫＦモデル及び定加速度ＫＦモデルそれぞれについて、次式（２４）で表される観測行列「Ｈ」を用いて計算を行う。

より具体的には、測位演算部３３は、図２のＩＭＭ演算処理のステップＳ１５において、観測行列「Ｈ」を用いてメジャメント予測値「ＨＸ⁻：Ｈ₁Ｘ⁻ ₁，Ｈ₂Ｘ⁻ ₂，Ｈ₃Ｘ⁻ ₃」を算出する。また、ステップＳ１７の補正演算において、観測行列「Ｈ」を用いて式（８）〜式（１０）に従って状態ベクトル及び誤差共分散の補正値「Ｘ：Ｘ₁，Ｘ₂，Ｘ₃」及び「Ｐ：Ｐ₁，Ｐ₂，Ｐ₃」を算出する。そして、ステップＳ１９及びＳ２１において、観測行列「Ｈ」を用いて式（１１）及び式（１２）に従って尤度指標値「Ｖｘ：Ｖｘ₁，Ｖｘ₂，Ｖｘ₃」及び重み係数「β：β₁，β₂，β₃」を算出する。

ここで、観測行列「Ｈ」に含まれる「Ｄ_x」、「Ｄ_y」及び「Ｄ_z」は、次式（２５−１）〜（２５−３）でそれぞれ与えられる。

但し、「ｘ_u」、「ｙ_u」及び「ｚ_u」は、それぞれＫＦ処理の予測演算により予測された携帯型電話機１の予測位置の３次元の成分をそれぞれ示しており、「ｘ_s」、「ｙ_s」及び「ｚ_s」は、観測対象とするＧＰＳ衛星の位置の３次元の成分をそれぞれ示している。また、「Ｒ」は、携帯型電話機１の予測位置と観測対象とするＧＰＳ衛星の位置間の距離を示している。

そして、測位演算部３３は、定点ＫＦモデル、等速度ＫＦモデル及び定加速度ＫＦモデルそれぞれについてステップＳ１９で算出した尤度指標値「Ｖｘ」の相対値「Ｖｘ₁−Ｖｘ₂」、「Ｖｘ₂−Ｖｘ₁」、「Ｖｘ₁−Ｖｘ₃」、「Ｖｘ₃−Ｖｘ₁」、「Ｖｘ₂−Ｖｘ₃」及び「Ｖｘ₃−Ｖｘ₂」を算出し（ステップＳ２５）、これらの相対値と、ステップＳ２１で算出した重み係数「β：β₁，β₂，β₃」と、ステップＳ３で算出した規格化定数「ＣＮ：Ｃ₁，Ｃ₂，Ｃ₃」とを用いて式（１３−２）に従ってモデル確率「μ：μ₁，μ₂，μ₃」を算出する（ステップＳ２７、Ｓ２９）。

そして、測位演算部３３は、算出したモデル確率「μ：μ₁，μ₂，μ₃」を用いて、式（１６）及び式（１７）に従って状態ベクトル「Ｘ：Ｘ₁，Ｘ₂，Ｘ₃」及び誤差共分散「Ｐ：Ｐ₁，Ｐ₂，Ｐ₃」を算出することで、状態ベクトルの出力値「Ｘ」及び誤差共分散の出力値「Ｐ」を決定する（ステップＳ３１）。

ステップＡ３でＩＭＭ演算処理を行った後、測位演算部３３は、ＩＭＭ演算処理で決定した状態ベクトルの出力値「Ｘ」に含まれる携帯型電話機１の３次元の位置「（ｘ，ｙ，ｚ）」を出力位置に決定し、ホストＣＰＵ５０に出力する（ステップＡ５）。

その後、測位演算部３３は、操作部６０を介してユーザにより測位終了指示がなされたか否かを判定し（ステップＡ７）、なされなかったと判定した場合は（ステップＡ７；Ｎｏ）、ステップＡ３に戻る。また、測位終了指示がなされたと判定した場合は（ステップＡ７；Ｙｅｓ）、ベースバンド処理を終了する。

４．作用効果
本実施形態によれば、カルマンフィルタ処理それぞれの尤度指標値「Ｖｘ」が算出され、この尤度指標値「Ｖｘ」のカルマンフィルタ処理間の相対値が算出される。そして、算出された相対値を用いてカルマンフィルタ処理それぞれのモデル確率「μ」が算出され、このモデル確率「μ」を用いてカルマンフィルタ処理それぞれの処理結果が合成されることで測位位置が決定される。

ここで、尤度指標値「Ｖｘ」とは、カルマンフィルタの入力値となる外部観測量の尤度「λ」を算出するための尤度関数の指数部分のことを指す。ＩＭＭ演算におけるモデル確率「μ」は、尤度指標値「Ｖｘ」を用いて算出した外部観測量の尤度「λ」を基に算出することができるが、尤度指標値「Ｖｘ」が非常に小さな値である場合には、計算機で計算しようとしても尤度「λ」の返し値が「０」になってしまい、モデル確率「μ」を正しく求めることができない。

この問題を解決するため、尤度指標値「Ｖｘ」をそのまま用いるのではなく、複数のカルマンフィルタ処理間の尤度指標値「Ｖｘ」の相対値を用いる。これは、個々の尤度指標値「Ｖｘ」が非常に小さな値であったとしても、尤度指標値「Ｖｘ」の相対値は概ね「０」の近傍の値となり、指数関数に当該相対値を与えて計算すると「０」とはならない点に着目したものである。かかる原理により、モデル確率「μ」を正しく算出することが可能となり、その結果、正確な測位位置の算出を実現し得る。

５．変形例
５−１．電子機器
本発明は、測位装置を備えた電子機器であれば何れの電子機器にも適用可能である。例えば、ノート型パソコンやＰＤＡ（Personal Digital Assistant）、カーナビゲーション装置等についても同様に適用可能である。

５−２．衛星測位システム
上述した実施形態では、衛星測位システムとしてＧＰＳを例に挙げて説明したが、ＷＡＡＳ（Wide Area Augmentation System）、ＱＺＳＳ（Quasi Zenith Satellite System）、ＧＬＯＮＡＳＳ（GLObal NAvigation Satellite System）、ＧＡＬＩＬＥＯ等の他の衛星測位システムであってもよい。

５−３．処理の分化
測位演算部３３が実行する処理の一部又は全部を、ホストＣＰＵ５０が実行することにしてもよい。例えば、メジャメント取得演算部３２から取得したメジャメント実測値を用いて、ホストＣＰＵ５０がＩＭＭ演算を行って出力位置を決定し、当該出力位置を表示部７０に表示させる。

ＩＭＭの説明図。 ＩＭＭ演算処理の流れを示すフローチャート。 ＩＭＭ演算処理の流れを示すフローチャート。 従来の手法によりモデル確率を求めた実験結果を示す図。 実施形態の手法によりモデル確率を求めた実験結果を示す図。 従来の手法によるモデル確率算出の説明図。 実施形態の手法によるモデル確率算出の説明図。 携帯型電話機の機能構成を示すブロック図。 ベースバンド処理回路部のＲＯＭに格納されたデータの一例を示す図。 ベースバンド処理回路部のＲＡＭに格納されるデータの一例を示す図。 実施例におけるＩＭＭの説明図。 捕捉衛星別メジャメントデータのデータ構成の一例を示す図。 ベースバンド処理の流れを示すフローチャート。

符号の説明

１ 携帯型電話機 、 １０ ＧＰＳアンテナ、 ２０ ＧＰＳ受信部、
２１ ＲＦ受信回路部、 ３０ ベースバンド処理回路部、 ３１ 演算制御部、
３２ メジャメント取得演算部、 ３３ 測位演算部、 ３５ ＲＯＭ、
３７ ＲＡＭ、 ４０ ＴＣＸＯ、 ５０ ホストＣＰＵ、 ６０ 操作部、
７０ 表示部、 ８０ 携帯電話用アンテナ、 ９０ 携帯電話用無線通信回路部、
１００ ＲＯＭ、 １１０ ＲＡＭ

Claims (6)

1. 測位装置が、各カルマンフィルタ処理（１，２，…ｊ）それぞれをモデルとみなし、前記カルマンフィルタ処理それぞれの処理結果を所定のモデル確率で重みづけして合成するインタラクティブミキシングモデル（Interactive Mixing Model）演算（以下「ＩＭＭ演算」と称す。）を行って測位する測位方法であって、
   当該カルマンフィルタ処理の状態ベクトル及び誤差共分散と所与のメジャメント実測値とに基づいて、前記カルマンフィルタ処理それぞれについての、当該メジャメント実測値の尤度を示す尤度関数の指数部分を示す尤度指標値（Ｖｘ _１ ，Ｖｘ _２ ，…Ｖｘ _ｊ ）を算出することと、
   前記算出された尤度指標値の前記カルマンフィルタ処理間の相対値（当該カルマンフィルタ処理ｊに関する相対値は（Ｖｘ _１ −Ｖｘ _ｊ ），（Ｖｘ _２ −Ｖｘ _ｊ ），…（Ｖｘ _ｎ −Ｖｘ _ｊ ）：但しｎはカルマンフィルタ処理の総数）を算出することと、
   以下の（式１）を用いて前記カルマンフィルタ処理それぞれの前記モデル確率（μ _１ ，μ _２ ，…μ _ｊ ）を算出することと、
   前記算出されたモデル確率を用いて前記カルマンフィルタ処理それぞれの処理結果を合成することで測位位置を決定することと、
   を含む測位方法。
   

   

   （但し、各カルマンフィルタ処理それぞれを１〜ｊとし、当該カルマンフィルタ処理の総数をｎとして添字で示す。Ｖｘが各モデルの尤度指標値、βが各モデルの重み係数、Ｃが各モデルの規格化定数を示す。）
2. 前記モデル確率を算出することは、前記相対値が所定の範囲内であるか否かを判定することを含み、前記所定の範囲内であると判定された場合に、前記（式１）を用いて前記モデル確率を算出することである、
   請求項１に記載の測位方法。
3. 前記モデル確率を算出することは、当該相対値が前記所定の範囲の上限以上のときには当該モデル確率を「１」、前記所定の範囲の下限以下のときには当該モデル確率を「１」とすることを含む、
   請求項２に記載の測位方法。
4. 前記複数のカルマンフィルタ処理には、前記測位装置の移動状態として、停止状態にあることを前提としたカルマンフィルタ処理と、等速移動状態にあることを前提としたカルマンフィルタ処理と、定加速度状態にあることを前提としたカルマンフィルタ処理との３種類のカルマンフィルタ処理が少なくとも含まれる、
   請求項１〜３の何れか一項に記載の測位方法。
5. 請求項１〜４の何れか一項に記載の測位方法を前記測位装置に内蔵されたコンピュータに実行させるためのプログラム。
6. 各カルマンフィルタ処理（１，２，…ｊ）それぞれをモデルとみなし、前記カルマンフィルタ処理それぞれの処理結果を所定のモデル確率で重みづけして合成するＩＭＭ演算を行って測位する測位装置であって、
   当該カルマンフィルタ処理の状態ベクトル及び誤差共分散と所与のメジャメント実測値とに基づいて、前記カルマンフィルタ処理それぞれについての、当該メジャメント実測値の尤度を示す尤度関数の指数部分を示す尤度指標値（Ｖｘ _１ ，Ｖｘ _２ ，…Ｖｘ _ｊ ）を算出する指標値算出部と、
   前記算出された尤度指標値の前記カルマンフィルタ処理間の相対値（（Ｖｘ _１ −Ｖｘ _ｊ ），（Ｖｘ _２ −Ｖｘ _ｊ ），…（Ｖｘ _ｎ −Ｖｘ _ｊ ）：但しｎはカルマンフィルタ処理の総数）を算出する相対値算出部と、
   以下の（式２）を用いて前記カルマンフィルタ処理それぞれの前記モデル確率（μ _１ ，μ _２ ，…μ _ｊ ）を算出するモデル確率算出部と、
   前記算出されたモデル確率を用いて前記カルマンフィルタ処理それぞれの処理結果を合成することで測位位置を決定する測位位置決定部と、
   を備えた測位装置。
   

   

   （但し、各カルマンフィルタ処理それぞれを１〜ｊとし、当該カルマンフィルタ処理の総数をｎとして添字で示す。Ｖｘが各モデルの尤度指標値、βが各モデルの重み係数、Ｃが各モデルの規格化定数を示す。）

Images