位置算出方法及び位置算出装置
【課題】位置算出の正確性を向上させる、位置算出方法及び位置算出装置を提供する。
【解決手段】測位用信号の一種であるＧＰＳ衛星信号を受信する。そして、少なくともＧＰＳ受信機の位置の変動を確率変数とする確率分布モデルである状態変動正規分布モデルに基づいて定められた位置算出演算を行って、ＧＰＳ受信機の位置を算出する。
【解決手段】測位用信号の一種であるＧＰＳ衛星信号を受信する。そして、少なくともＧＰＳ受信機の位置の変動を確率変数とする確率分布モデルである状態変動正規分布モデルに基づいて定められた位置算出演算を行って、ＧＰＳ受信機の位置を算出する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、位置算出方法及び位置算出装置に関する。
【背景技術】
【０００２】
測位用信号を利用した測位システムとしては、ＧＰＳ（Global Positioning System）が広く知られており、携帯型電話機やカーナビゲーション装置等に内蔵された位置算出装置に利用されている。ＧＰＳでは、複数のＧＰＳ衛星の位置や各ＧＰＳ衛星から位置算出装置までの擬似距離等の情報に基づいて、位置算出装置の位置座標及び時計誤差を求める位置算出演算を行う。
【０００３】
衛星測位システムを利用した位置算出演算では、複数の測位用衛星について受信信号を用いて測定した擬似距離に含まれ得る誤差（以下、「擬似距離誤差」と称す。）の二乗和を最小化させる、いわゆる最小二乗法を利用した位置算出演算が知られている（例えば、特許文献１。）。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００９−９７８９７号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
従来の位置算出演算では、位置算出演算を行う時刻（以下、「演算時刻」と称す。）毎に、最小二乗法を利用するなどして位置算出演算を１回１回個別に行うのが通常である。従って、例えば、測定した擬似距離に瞬間的に大きな誤差が混入した場合、そのときの位置算出演算の精度が急激に低下するという欠点があったが、以降の位置算出演算には影響しないという利点もあった。但し、瞬間的とは言え、真の位置から大きく離れた位置が演算結果として得られる、いわゆる位置飛びが発生することは、場合によって大きな問題となり得た。
【０００６】
本発明は上述した課題に鑑みて為されたものであり、その目的とするところは、位置算出の正確性を向上させるための新しい位置算出手法を提案することにある。
【課題を解決するための手段】
【０００７】
以上の課題を解決するための第１の形態は、測位用信号を受信することと、前記測位用信号の受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行うことと、を含む位置算出方法である。
【０００８】
また、他の形態として、測位用信号を受信する受信部と、前記受信部による受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う位置算出部と、を備えた位置算出装置を構成してもよい。
【０００９】
この第１の形態等によれば、測位用信号の受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う。これにより、各演算結果を独立とするのではなく、演算結果に時間的な相関を持たせて位置算出を行うことが可能となる。この場合、測定量に瞬間的に大きな誤差が混入したとしても、位置算出精度が低下することを防止できる。
【００１０】
また、第２の形態として、第１の形態の位置算出方法において、過去の前記位置算出演算の結果に基づいて前記変動の基準を更新することと、前記基準を用いて前記変動を算出することと、を更に含む、位置算出方法を構成することとしてもよい。
【００１１】
この第２の形態によれば、過去の位置算出演算の結果に基づいて演算結果の変動の基準を更新することで、過去の位置算出演算の結果を今回の位置算出演算の結果に反映させることが可能となる。更新した基準を用いて演算結果の変動を算出することで、位置算出の正確性が向上する。
【００１２】
また、第３の形態として、第１又は第２の形態の位置算出方法において、前記位置算出演算は、少なくとも前記演算結果の変動を確率変数とした正規分布モデルに基づいて定められた演算であり、過去の前記位置算出演算の結果に基づいて前記正規分布モデルを更新することを更に含む、位置算出方法を構成することとしてもよい。
【００１３】
この第３の形態によれば、少なくとも演算結果の変動を確率変数とした正規分布モデルに基づいて定められた位置算出演算を行う。つまり、演算結果の変動の分布が正規分布に従うと仮定する。その上で、正規分布モデルを固定的に設定するのではなく、過去の位置算出演算の結果に基づいて正規分布モデルを更新することで、時間変化に対応した適切な正規分布モデルに基づく位置算出が可能となる。
【００１４】
また、第４の形態として、第３の形態の位置算出方法において、前記受信信号に基づいて擬似距離を測定することを更に含み、前記位置算出演算は、前記演算結果の変動を確率変数とした正規分布モデルと、前記擬似距離に含まれる誤差を確率変数とした確率分布モデルとに基づいて定められた演算である、位置算出方法を構成することとしてもよい。
【００１５】
この第４の形態によれば、測位用信号の受信信号に基づいて擬似距離を測定する。そして、演算結果の変動を確率変数とした正規分布モデルと、擬似距離に含まれる誤差を確率変数とした確率分布モデルとに基づいて定められた位置算出演算を行う。演算結果の変動を確率変数とした正規分布モデルに加えて、擬似距離に含まれる誤差を確率変数とした確率分布モデルを用いることで、位置算出演算を適切に行うことができる。
【００１６】
より具体的には、第５の形態として、第４の形態の位置算出方法における前記擬似距離に含まれる誤差を確率変数とした確率分布モデルを、正規分布モデル又は正規混合分布モデルとする、位置算出方法を構成することができる。
【００１７】
通常の測位環境においては、擬似距離に含まれる誤差の分布は正規分布に従うと仮定することができる。しかし、例えばマルチパス環境のように、擬似距離に含まれる誤差の分布を正規分布と仮定することが必ずしも適切でない場合もある。そこで、第５の形態のように、擬似距離に含まれる誤差を確率変数とした確率分布モデルを、正規分布モデル又は正規混合分布モデルとして位置算出演算を行うことで、測位環境に見合った適切な位置算出を実現することが可能となる。
【００１８】
また、第６の形態として、第３〜第５の何れかの形態の位置算出方法において、前記位置算出演算は、位置情報の他に、クロックバイアス、速度情報及びクロックドリフトのうちの何れかを算出する演算であり、前記正規分布モデルは、当該演算結果の各要素の変動を確率変数としたモデルである、位置算出方法を構成することとしてもよい。
【００１９】
この第６の形態によれば、位置情報の他に、クロックバイアス、速度情報及びクロックドリフトのうちの何れかを算出する位置算出演算を行う。その際、演算結果の各要素の変動を確率変数とした正規分布モデルを用いることで、各要素に時間的な相関を持たせて位置算出を行うことが可能となる。
【００２０】
また、第７の形態として、第１〜第６の何れかの形態の位置算出方法において、前記位置算出演算によって算出可能な算出パラメーターのうちの所定のパラメーターの値を、前記受信信号を用いた所定の収束演算を行って算出することを更に含み、前記位置算出演算は、過去に前記収束演算によって算出された前記所定のパラメーターの値を用いる演算であることを含む、位置算出方法を構成することとしてもよい。
【００２１】
この第７の形態によれば、位置算出演算によって算出可能な算出パラメーターのうちの所定のパラメーターの値を、受信信号を用いた所定の収束演算を行って算出する。そして、過去に収束演算によって算出された所定のパラメーターの値を用いる位置算出演算を行う。これにより、特定の算出パラメーターについては、時間的な相関を持たせずに、所定の収束演算を行って値を算出することができる。そして、その算出結果を、所与の確率分布モデルに基づいて定められた位置算出演算に反映させることで、位置算出を効果的に行うことが可能となる。
【図面の簡単な説明】
【００２２】
【図１】状態変動正規分布モデルの説明図。
【図２】擬似距離誤差正規分布モデルの説明図。
【図３】マルチパスの説明図。
【図４】自己相関の一例を示す図。
【図５】マルチパス環境における自己相関の一例を示す図。
【図６】マルチパス環境における自己相関の一例を示す図。
【図７】第１の擬似距離誤差正規混合分布モデルの説明図。
【図８】第２の擬似距離誤差正規混合分布モデルの説明図。
【図９】位置算出結果の一例を示す図。
【図１０】位置算出結果の一例を示す図。
【図１１】位置算出結果の一例を示す図。
【図１２】携帯型電話機の機能構成の一例を示すブロック図。
【図１３】ベースバンド処理回路部の回路構成の一例を示す図。
【図１４】ベースバンド処理の流れを示すフローチャート。
【図１５】第１の位置算出処理の流れを示すフローチャート。
【図１６】第２の位置算出処理の流れを示すフローチャート。
【図１７】第３の位置算出処理の流れを示すフローチャート。
【図１８】第４の位置算出処理の流れを示すフローチャート。
【図１９】ドップラー測位処理の流れを示すフローチャート。
【発明を実施するための形態】
【００２３】
以下、図面を参照して、本発明の好適な実施形態の一例について説明する。本実施形態は、衛星測位システムの一種であるＧＰＳ（Global Positioning System）を適用した実施形態である。なお、本発明を適用可能な実施形態が以下説明する実施形態に限定されるわけでないことは勿論である。
【００２４】
１．原理
本実施形態における位置算出方法（位置算出の原理）について説明する。ＧＰＳを利用した衛星測位システムにおいて、測位用衛星の一種であるＧＰＳ衛星は、エフェメリスやアルマナックといった衛星軌道データを含む航法データを、測位用の衛星信号の一種であるＧＰＳ衛星信号に乗せて発信している。
【００２５】
ＧＰＳ衛星信号は、拡散符号の一種であるＣ／Ａ（Coarse and Acquisition）コードによって、スペクトラム拡散方式として知られるＣＤＭＡ（Code Division Multiple Access）方式によって変調された１．５７５４２［ＧＨｚ］の通信信号である。Ｃ／Ａコードは、コード長１０２３チップを１ＰＮフレームとする繰返し周期１ｍｓの擬似ランダム雑音符号であり、各ＧＰＳ衛星に固有のコードである。
【００２６】
ＧＰＳでは、複数のＧＰＳ衛星の位置や各ＧＰＳ衛星からＧＰＳ受信機（位置算出装置）までの擬似距離といった、ＧＰＳ衛星信号の受信信号に係る測定量を用いて、ＧＰＳ受信機の位置（位置座標）や時計誤差（クロックバイアス）を演算する位置算出演算を行う。擬似距離は、衛星サーチを行ってＧＰＳ衛星信号を捕捉し、その結果として得られるメジャメント情報を用いて測定される。
【００２７】
衛星サーチには、ＧＰＳ衛星信号を受信した信号に対する位相方向の相関演算（以下、「位相サーチ」と称す。）と、周波数方向の相関演算（以下、「周波数サーチ」と称す。）とが含まれる。位相サーチを行うことで、ＧＰＳ受信機がＧＰＳ衛星信号を受信した際のＣ／Ａコードの位相（以下、「コード位相」と称す。）をメジャメント情報として取得できる。また、周波数サーチを行うことで、ＧＰＳ受信機が受信したＧＰＳ衛星信号の受信周波数やドップラー周波数をメジャメント情報として取得できる。
【００２８】
観念的には、ＧＰＳ衛星とＧＰＳ受信機との間には複数のＣ／Ａコードが並んでいると考えることができる。ＧＰＳ衛星とＧＰＳ受信機との間の距離は、ちょうどＣ／Ａコードの整数倍の長さになるとは限らず、端数部分が生じ得る。この擬似距離の端数部分に相当するのがコード位相である。擬似距離の整数部分は、ＧＰＳ受信機及びＧＰＳ衛星の概略位置から求められる。
【００２９】
本実施形態では、ＧＰＳ受信機（位置算出装置）の位置ベクトル“Ｐ”及びクロックバイアス（時計誤差）“ｂ”を未知数とし、これらを成分とするＧＰＳ受信機の状態ベクトル“Ｘ”を次式（１）のように定義する。
【数１】
【００３０】
式（１）において、“（ｘ，ｙ，ｚ）”は、測位を行う座標系におけるＧＰＳ受信機（位置算出装置）の位置ベクトル“Ｐ”の各軸の位置成分である。後々の便宜のため、位置成分“（ｘ，ｙ，ｚ）”及びクロックバイアス“ｂ”を、状態“Ｘ”の成分という意味で“Ｘi＝（ｘ，ｙ，ｚ，ｂ）T＝（Ｘ1，Ｘ2，Ｘ3，Ｘ4）T”と表記する。但し、下付きの添え字の“ｉ＝１，２，３，４”は、状態“Ｘ”の成分の番号を示す。上付きの添え字の“Ｔ”は行列の転置を示す。また、本実施形態の数式では、ベクトルを意味する矢印の表記を省略して簡潔に記載する。
【００３１】
ＧＰＳ受信機の真の位置ベクトルを“Ｐ＝Ｐ0＋δＰ”、真のクロックバイアスを“ｂ＝ｂ0＋δｂ”とおく。但し、“δＰ”及び“δｂ”は、それぞれ位置ベクトル及びクロックバイアスの初期値に適用される未知の変化量（修正量）であり、それぞれ「位置変化量ベクトル」及び「クロックバイアス変化量」と呼称する。また、位置変化量ベクトル“δＰ”及びクロックバイアス変化量“δｂ”を纏めて、次式（２）のように状態変化量ベクトル“δＸ”を定義する。
【数２】
【００３２】
従来から用いられているＧＰＳの測位計算では、例えば次式（３）の観測方程式を利用して、状態変化量ベクトル“δＸ”を求める。
【数３】
【００３３】
本明細書では、ＧＰＳ衛星のうちＧＰＳ受信機が捕捉に成功した衛星（以下、「捕捉衛星」と称す。）の番号を“ｋ”と表記し、原則的に各諸量について下付きの添え字で表記する。また、本実施形態では、ＧＰＳ受信機が合計Ｋ個のＧＰＳ衛星の捕捉に成功したものとして説明する。つまり、「ｋ＝１，２，・・・，Ｋ」として説明する。
【００３４】
式（３）の左辺の“δρ”は擬似距離修正量ベクトルであり、次式（４）のように表される。
【数４】
つまり、擬似距離修正量ベクトル“δρ”は、各捕捉衛星の擬似距離修正量“δρk”を成分とするベクトルである。
【００３５】
また、式（３）の右辺の行列“Ｇ”は、捕捉衛星の衛星配置を決定付ける幾何行列であり、次式（５）ように表される。
【数５】
但し、“ｌk”は、ＧＰＳ受信機からｋ番目の捕捉衛星に向かう視線方向のベクトル（以下、「視線ベクトル」と称す。）である。
【００３６】
また、式（３）の右辺の“Ｅ”は擬似距離誤差ベクトルであり、次式（６）のように表される。
【数６】
但し、“εk”は、ｋ番目の捕捉衛星からＧＰＳ受信機までの擬似距離に含まれる誤差（以下、「擬似距離誤差」と称す。）である。擬似距離誤差ベクトル“Ｅ”は、各捕捉衛星の擬似距離誤差“εk”を成分とするベクトルである。
【００３７】
式（３）における未知数は、位置変化量ベクトル“δＰ＝（δｘ，δｙ，δｚ）”及びクロックバイアス変化量“δｂ”の合計４個である。従って、式（３）は、Ｋ≧４の場合に解くことができる。Ｋ＞４の場合は、式（３）は過決定方程式となる。この場合は、Ｋ個の捕捉衛星の擬似距離誤差“ε”の二乗和を最小とする解を求める最小二乗法が一般的に用いられる。具体的には、式（３）の最小二乗解は、例えば次式（７）のように求められる。
【数７】
【００３８】
式（３）の観測方程式を利用した位置算出演算では、演算時刻毎に、ＧＰＳ衛星信号を受信して擬似距離を測定し、式（７）の最小二乗解を解析的に求める。この場合、状態変化量ベクトル“δＸ”を、演算時刻毎に独立として計算している。
【００３９】
しかし、ＧＰＳ受信機（位置算出装置）は、通常は連続的に移動する。例えば、ＧＰＳ受信機を携行したユーザーが通常歩行や通常走行しているのであれば、演算時刻毎のＧＰＳ受信機の移動距離にそれほど大きな差は現れないはずである。つまり、ユーザーが突飛な動作をしない限り、ＧＰＳ受信機の位置変化量ベクトル“δＰ”の変動はほぼ一定となることが期待される。ＧＰＳ受信機を自動車等の移動体に設置した場合も同様である。このことから、本願発明者は、位置算出演算の演算結果の変動を確率変数とした所与の確率分布モデルを定義することを考えた。
【００４０】
例えば、ユーザーが通常歩行しており、「１秒間隔」で位置算出演算を行うとすると、１秒毎のＧＰＳ受信機の移動距離は「数メートル程度」となり、個人差を考慮しても、その誤差幅は「数十センチ程度」になると考えられる。従って、歩行を想定した場合、位置算出演算の演算結果のうち、位置変化量ベクトル“δＰ”については、期待値が「数メートル程度」、標準偏差が「数十センチ程度」と仮定できるであろう。また、ＧＰＳ受信機を自動車に設置する場合を考えると、位置変化量ベクトル“δＰ”は、期待値が「数十メートル程度」、標準偏差が「数メートル程度」と仮定できるであろう。
【００４１】
クロックバイアス“ｂ”についても同様に考えると、短い時間間隔で考えた場合、クロックバイアス“ｂ”がそれほど大きく変動するとは考えにくい。そのため、クロックバイアス変化量“δｂ”についても確率分布モデルを仮定できるであろう。
【００４２】
そこで、本願発明者は、位置ベクトル“Ｐ”及びクロックバイアス“ｂ”を成分とする状態ベクトル“Ｘ”について、その変動である状態変化量ベクトル“δＸ”を確率変数とする正規分布型の確率分布モデル（以下、「状態変動正規分布モデル」と称す。）を仮定した。
【００４３】
図１は、状態変動正規分布モデル（状態変動ＧＭ（Gaussian Model））の説明図である。図１には、状態変化量ベクトル“δＸ”の１つの成分“δＸi”に着目したモデルを図示している。横軸は状態変化量の成分“δＸi”であり、縦軸はその確率密度“ｐ（δＸi）”である。状態変化量“δＸi”の期待値“μ”を「状態変化量期待値」と呼称し、状態変化量“δＸi”の標準偏差“σ”を「状態変化量標準偏差」と呼称する。
【００４４】
状態変化量“δＸi”を確率変数とし、その分布が正規分布Ｎ（μ，σ2）に従うと仮定する。しかし、本実施形態では、状態変化量ベクトル“δＸ”には、位置変化量ベクトル“Ｐ”の各成分“（δｘ，δｙ，δｚ）”及びクロックバイアス変化量“δｂ”が含まれる。そこで、一次元の正規分布モデルを拡張して、次式（８）で与えられる多次元状態変動正規分布モデル“ｐδX(m)（δＸ）”を用いることとする。
【数８】
【００４５】
式（８）において、上付きの括弧書きの“ｍ”は確率変数として取り扱う算出パラメーターの数を示す。算出パラメーターとは、位置算出演算によって算出可能なパラメーターであり、本実施形態では“（δｘ，δｙ，δｚ，δｂ）”の４つのパラメーターがこれに相当する。４つの算出パラメーター全てを確率変数とするのであれば“ｍ＝４”となる。詳細は後述するが、４つの算出パラメーターのうちの一部の算出パラメーターを確率変数とすることも可能である。この場合は、確率変数とする算出パラメーターの数が“ｍ”となる。
【００４６】
式（８）の多次元状態変動正規分布モデル“ｐδX(m)（δＸ）”は、非独立のｍ個の確率変数の同時発生確率を与える確率密度関数である。“μδX”は、状態変化量ベクトル“δＸ”の各成分の期待値でなる状態変化量期待値ベクトルである。また、“ＳδX”は、状態変化量ベクトル“δＸ”の成分間の共分散でなる状態変化量共分散行列である。これらについては、数式を用いて詳細に後述する。
【００４７】
本実施形態では、上記の多次元状態変動正規分布モデル“ｐδX(m)（δＸ）”を適用して、最尤推定法を利用して状態変化量ベクトル“δＸ”の最尤推定解を求める。これは、測位用信号を受信し、測位用信号の受信信号を用いた演算であって、位置算出演算の演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行うことに相当する。状態変化量ベクトル“δＸ”の最尤推定解が求まれば、当該最尤推定解を状態ベクトルの初期値“Ｘ0”に加算することで、真の状態ベクトル“Ｘ”が推定できる。
【００４８】
さて、本実施形態の位置算出演算では、ＧＰＳ衛星信号の受信信号に基づいて擬似距離を測定する。そして、上記の状態変動正規分布モデルと、擬似距離に含まれる誤差を確率変数とした確率分布モデルとに基づいて定められた位置算出演算を行う。つまり、擬似距離誤差“ε”も確率変数として取り扱い、所与の確率分布モデルを仮定して演算を行う。本実施形態では、擬似距離誤差“ε”を確率変数とした確率分布モデルを、正規分布モデル又は正規混合分布モデルとする。
【００４９】
図２は、擬似距離誤差正規分布モデル（擬似距離誤差ＧＭ）の説明図である。図２において、横軸は擬似距離誤差“ε”であり、縦軸はその確率密度“ｐε（ε）”である。擬似距離誤差“ε”の期待値“με”を「擬似距離誤差期待値」と呼称し、擬似距離誤差“ε”の標準偏差“σε”を「擬似距離誤差標準偏差」と呼称する。
【００５０】
擬似距離誤差“ε”が、正規分布Ｎ（με，（σε）2）に従うと仮定して観測方程式の解を求める手法が最小二乗法である。広く用いられるのは「με＝０，σε＝１」とする擬似距離誤差正規分布モデルであり、この場合の最小二乗解が式（７）と一致する。
【００５１】
上記の擬似距離誤差正規分布モデル“ｐε（ε）”は、ＧＰＳ受信機の天空が開けた環境（例えばオープンスカイ環境）において有用である。オープンスカイ環境では、ＧＰＳ衛星から品質の良いＧＰＳ衛星信号を受信することができることから、擬似距離誤差“ε”はゼロ近傍の値となることが期待されるためである。しかし、ＧＰＳ受信機の天空が開けていない環境では、擬似距離誤差“ε”の挙動が変化する。その典型例はマルチパス環境である。
【００５２】
図３は、マルチパス環境の説明図である。直接波信号を点線で示し、間接波信号を一点鎖線で示す。間接波信号の存在により、ＧＰＳ受信機が測定する擬似距離には誤差が生ずる。擬似距離は、ＧＰＳ受信機が位相サーチを行って取得したコード位相を用いて算出されるが、マルチパス環境では、このコード位相に誤差が重畳する。コード位相の誤差は、真のコード位相に対して、正の誤差となる場合及び負の誤差となる場合がある。
【００５３】
図４は、コード位相検出の原理の説明図である。図４では、横軸をコード位相、縦軸を相関値として、Ｃ／Ａコードの自己相関値の概略例を示している。なお、以下の説明では、相関値というときは、相関値の大きさ（絶対値）を意味するものとする。
【００５４】
Ｃ／Ａコードの自己相関値は、理想的にはピーク値を頂点とする左右対称の略三角形の形状で表される。この場合に、相関値のピーク値（以下、「相関ピーク値」と称す。）に対応する位相が、受信したＧＰＳ衛星信号のＣ／Ａコードの位相である。相関ピーク値を検出するために、例えば、あるコード位相に対して、一定量だけ進んだコード位相（以下、「進み位相」と称す。）と、一定量だけ遅れたコード位相（以下、「遅れ位相」と称す。）とを考える。そして、進み位相の相関値と遅れ位相の相関値とが等しくなる相関値を相関ピーク値として検出する。
【００５５】
図４では、相関値は左右対称の略三角形の形状であるため、進み位相の相関値と遅れ位相の相関値との中心のコード位相の相関値が相関ピーク値となり、対応するコード位相が「ピーク位相」として検出される。この検出されたピーク位相のことを「検出ピーク位相」と称する。図４の相関値の形状は理想形状であるが、マルチパス環境では相関値の形状が変化する。
【００５６】
図５及び図６は、マルチパス環境における相関値の形状の一例を示す図である。図５は、間接波信号が直接波信号と同位相で到達した場合の相関値のグラフの一例であり、図６は、間接波信号が直接波信号と逆位相で到達した場合の相関値のグラフの一例である。これらの図には、直接波信号と、間接波信号と、この直接波信号と間接波信号とを合成したマルチパス信号とのそれぞれに対応する相関値のグラフを示している。横軸はコード位相、縦軸は相関値である。
【００５７】
間接波信号に対する相関値は、直接波信号に対する相関値と同様に略三角形の形状をなしているが、間接波信号の相関ピーク値の大きさは、直接波信号の相関ピーク値よりも小さい。これは、ＧＰＳ衛星から送出されたＧＰＳ衛星信号が、建物や地面に反射したり障害物を透過することによって、送出された時点における信号強度が、受信時には弱められていることによるものである。
【００５８】
また、間接波信号のピーク位相は、直接波信号のピーク位相よりも遅れている。これは、ＧＰＳ衛星から送出されたＧＰＳ衛星信号が、建物や地面に反射したり障害物を回折することによって、ＧＰＳ衛星からＧＰＳ受信機までの伝搬距離が長くなったことによるものである。マルチパス信号に対する相関値は、直接波信号の相関値と間接波信号の相関値との和となるため、三角形状が歪んでピーク値を中心とした左右対称とはならない。
【００５９】
間接波信号が直接波信号と同位相でＧＰＳ受信機に到達した場合は、直接波信号と間接波信号とは互いに強め合う。このため、合成波信号の相関値は、直接波信号に対する相関値の大きさと間接波信号に対する相関値の大きさとの合算値となる。この場合、相関値の形状は、例えば図５に示すような形状となる。図５において、検出ピーク位相は、真のピーク位相よりも遅れた位相となる。
【００６０】
それに対し、間接波信号が直接波信号よりも、例えば半周期以上１周期未満の範囲等の逆位相となった場合は、直接波信号と間接波信号とは互いに弱め合う。このため、合成波信号の相関値は、直接波信号に対する相関値の大きさから間接波信号に対する相関値の大きさを減じた値となる。この場合、相関値の形状は、例えば図６に示すような形状となる。図６において、検出ピーク位相は、真のピーク位相よりも進んだ位相となる。
【００６１】
ここで、検出ピーク位相と真のピーク位相との位相差を「コード位相誤差」と定義する。そして、便宜的に、検出ピーク位相が真のピーク位相よりも遅れている場合のコード位相誤差の符号を「正」、検出ピーク位相が真のピーク位相よりも進んでいる場合のコード位相誤差の符号を「負」と定義する。この場合、図５ではコード位相誤差は「正」となり、図６ではコード位相誤差は「負」となる。
【００６２】
前述したように、擬似距離はコード位相を用いて算出されるため、コード位相の誤差は擬似距離に誤差として重畳される。すなわち、図５のようにコード位相に正の誤差が含まれると、擬似距離誤差“ε”は正の誤差となる。一方、図６のようにコード位相に負の誤差が含まれると、擬似距離誤差“ε”は負の誤差となる。
【００６３】
このように、マルチパス環境等においては、コード位相誤差の正負の符号に応じて、擬似距離誤差“ε”は正の値にもなり得るし、負の値にもなり得る。勿論、擬似距離誤差“ε”にはコード位相誤差以外の誤差成分も含まれる。そのため、これらの誤差成分の大きさ如何によっても、擬似距離誤差“ε”は種々の値をとり得る。
【００６４】
このことから、本願発明者は、例えばマルチパス環境においては、擬似距離誤差“ε”の分布を正規分布と仮定するのではなく、複数の正規分布を混合した正規混合分布と仮定することが適切であると判断した。擬似距離誤差“ε”に正規混合分布を仮定した確率分布モデルとして、本実施形態では２種類のモデルを例示する。
【００６５】
図７は、擬似距離誤差ＧＭＭの一例である第１の擬似距離誤差正規混合分布モデル（第１の擬似距離誤差ＧＭＭ（Gaussian Mixture Model））の説明図である。第１の擬似距離誤差ＧＭＭ“ｆ（ε）”は、第１の正規分布モデル“ｐ1（ε）”と、第２の正規分布モデル“ｐ2（ε）”との２種類の正規分布モデルを合成した関数として定義される。図７において、横軸は擬似距離誤差“ε”（単位はメートル）であり、縦軸はその確率密度“ｆ（ε）”である。また、第１の正規分布モデル“ｐ1（ε）”を点線で示し、第２の正規分布モデル“ｐ2（ε）”を一点鎖線で示し、第１の擬似距離誤差ＧＭＭ“ｆ（ε）”を太実線でそれぞれ示す。
【００６６】
第１の正規分布モデル“ｐ1（ε）”は、ＧＰＳ受信機が直接波信号を受信することを想定した擬似距離誤差“ε”の確率分布モデルである。間接波信号の影響がなければ、擬似距離誤差“ε”は“０メートル”の近傍の値となることが期待される。そのため、擬似距離誤差期待値を“０メートル”（με＝０）とし、擬似距離誤差標準偏差を比較的小さな値である“１０メートル”（σε＝１０）とする正規分布関数“Ｎ（０，１０2）”によって、第１の正規分布モデル“ｐ1（ε）”をモデル化した。
【００６７】
第２の正規分布モデル“ｐ2（ε）”は、ＧＰＳ受信機が直接波信号及び間接波信号を受信することを想定した擬似距離誤差“ε”のモデル関数である。前述したように、コード位相誤差の正負の符号により、擬似距離誤差“ε”は正負の何れの値もとり得る。そのため、第２の正規分布モデル“ｐ2（ε）”の擬似距離誤差期待値を、正負の中心である“０メートル”（με＝０）とした。また、コード位相誤差の大きさ如何によっては、擬似距離誤差“ε”に非常に大きな誤差が含まれ得る。そのため、擬似距離誤差標準偏差を大きめに見積もって“１００メートル”（σε＝１００）とした。すなわち、第２の正規分布モデル“ｐ2（ε）”を正規分布関数“Ｎ（０，１００2）”によってモデル化した。
【００６８】
第１の正規分布モデル“ｐ1（ε）”と第２の正規分布モデル“ｐ2（ε）”とを混合する割合（合成比率）は、次のように決定した。ＧＰＳ受信機がマルチパスの影響を受ける割合は、一般的にはそれほど高くないと考えられる。そこで、ＧＰＳ受信機が１割の確率でマルチパスの影響を受けると仮定して、第１の正規分布モデル“ｐ1（ε）”の重みを“９／１０”（９割）とし、第２の正規分布モデル“ｐ1（ε）”の重みを“１／１０”（１割）とした。従って、第１の正規分布モデル“ｐ1（ε）”を９／１０倍した関数と、第２の正規分布モデル“ｐ2（ε）”を１／１０倍した関数とを混合（合成）した関数として、第１の擬似距離誤差ＧＭＭ“ｆ（ε）”をモデル化した。
【００６９】
図８は、擬似距離誤差ＧＭＭの別例である第２の擬似距離誤差正規混合分布モデル（第２の擬似距離誤差ＧＭＭ）の説明図である。第２の擬似距離誤差ＧＭＭ“ｆ（ε）”は、図７の第１の擬似距離誤差ＧＭＭと同様に、第１の正規分布モデル“ｐ1（ε）”と第２の正規分布モデル“ｐ2（ε）”との２種類の正規分布モデルを混合（合成）した確率密度関数である。
【００７０】
第１の正規分布モデル“ｐ1（ε）”は、ＧＰＳ受信機が直接波信号を受信することを想定した擬似距離誤差“ε”のモデル関数として、正規分布関数“Ｎ（０，１０2）”でモデル化した。また、第２の正規分布モデル“ｐ2（ε）”は、マルチパス環境においてＧＰＳ受信機が間接波信号を受信することを想定した擬似距離誤差“ε”のモデル関数として、正規分布関数“Ｎ（１００，２０2）”でモデル化した。
【００７１】
図７の第１の擬似距離誤差ＧＭＭ“ｆ（ε）”と異なるのは、第２の正規分布モデル“ｐ2（ε）”を、直接波信号及び間接波信号を受信することを想定したモデル関数として定義するのではなく、間接波信号単体を受信することを想定したモデル関数として定義した点にある。これは、例えば、ＧＰＳ受信機の近傍に高層ビル等の障害物が存在する環境において起こり得る。ＧＰＳ受信機の天空が完全に開けていないために、ＧＰＳ衛星から送出された信号は障害物に反射し、間接波信号としてＧＰＳ受信機に到達する。しかし、障害物に遮られて、直接波信号がＧＰＳ受信機に直接到達しない環境である。アーバンキャニオン環境が典型例である。
【００７２】
本願発明者が実験を行った結果、アーバンキャニオン環境では、擬似距離誤差“ε”として「＋１００〜＋１５０メートル」程度の正の誤差が含まれる傾向があることがわかった。そこで、第２の擬似距離誤差ＧＭＭ“ｆ（ε）”では、第２の正規分布モデル関数“ｐ2（ε）”の擬似距離誤差期待値を“１００メートル”（με＝１００）とし、擬似距離誤差標準偏差を“２０メートル”（σε＝２０）とした。すなわち、“ｐ2（ε）＝Ｎ（１００，２０2）”としてモデル化した。
【００７３】
第１及び第２の正規分布モデルの合成比率は、例えばＧＰＳ受信機が２割の確率でマルチパスの影響を受けると仮定して決定した。すなわち、第１の正規分布モデル“ｐ1（ε）”の重みを“８／１０”（８割）とし、第２の正規分布モデル“ｐ2（ε）”の重みを“２／１０”（２割）とした。従って、第２の正規分布モデル“ｐ1（ε）”を８／１０倍した関数と、第２の正規分布モデル“ｐ2（ε）”を２／１０倍した関数とを混合（合成）した関数として、第２の擬似距離誤差ＧＭＭ“ｆ（ε）”をモデル化した。
【００７４】
上記のように、擬似距離誤差“ε”を確率変数とする確率分布モデルとして、正規分布モデル及び正規混合分布モデルの２種類のモデルを定義した。位置算出を行う際には、測位環境に応じて、擬似距離誤差正規分布モデルと、擬似距離誤差正規混合分布モデルとを使い分けて位置算出演算を行う。以下、それぞれの確率分布モデルを適用する場合の具体的な位置算出方法について、数式を用いて詳細に説明する。
【００７５】
１−１．擬似距離誤差正規分布モデル（擬似距離誤差ＧＭ）を適用する場合
（１）第１の位置算出方法
第１の位置算出方法は、擬似距離誤差正規分布モデルを適用し、状態変化量ベクトル“δＸ”の全ての算出パラメーターを確率変数として位置算出する方法である。この方法では、状態変化量ベクトル“δＸ”の４つの算出パラメーター“（δｘ，δｙ，δｚ，δｂ）”全てを確率変数とする。この場合、式（８）の多次元状態変動正規分布モデル“ｐδX(m)（δＸ）”は、４次元状態変動正規分布モデル“ｐδX(4)（δＸ）”となる。
【００７６】
最尤推定法は、尤度と呼ばれる指標値に基づいて、所与の観測データを利用して推定対象とするパラメーターの最尤解を推定する統計学的手法である。尤度は、ある観測データについて、適用する確率モデルがどの程度当てはまっているかを確率で表す尺度である。
【００７７】
本実施形態において、尤度Ｌは、式（６）の擬似距離誤差ベクトル“Ｅ”のＫ衛星の擬似距離誤差“εk”（ｋ＝１，２，・・・，Ｋ）の同時発生確率と、状態変化量ベクトル“δＸ”の４つの算出パラメーター“δＸi＝（δｘ，δｙ，δｚ，δｂ）”（ｉ＝１，２，３，４）の同時発生確率との積で表される。つまり、尤度Ｌは、次式（９）のように表される。
【数９】
【００７８】
式（９）において、“ｐε,k（ε）”はｋ番目の捕捉衛星の擬似距離誤差正規分布モデルを示す。なお、図２で定義した擬似距離誤差正規分布モデル“ｐε（ε）”は、全ての捕捉衛星についてパラメーター値を共通化した共通のモデルとしてもよいし、捕捉衛星毎にパラメーター値を個別化した異なるモデルとしてもよい。
【００７９】
また、式（９）において、“ｐδX,i（δＸi）”はｉ番目の算出パラメーターの状態変動正規分布モデルを示す。１行目から２行目に変形する際に、４つの算出パラメーターの同時発生確率を与える４次元状態変動正規分布モデル“ｐδX(4)（δＸ）”に書き換えた。
【００８０】
式（９）の尤度Ｌを最大化する状態変化量ベクトル“δＸ”を最尤推定解として求める。尤度Ｌをそのまま用いて計算するのは困難であるため、本実施形態では、その対数をとった対数尤度ｌｏｇＬを用いて計算する。対数尤度ｌｏｇＬは、次式（１０）のように計算される。
【数１０】
【００８１】
式（１０）において、“ｃｏｎｓｔ”は定数項である。また、“μE”は、Ｋ衛星の擬似距離誤差“ε”の期待値“με,k＝（με,1，με,2，・・・，με,K）”を成分とする擬似距離誤差期待値ベクトルであり、次式（１１）で与えられる。
【数１１】
【００８２】
“ＷE”は、擬似距離誤差逆共分散行列であり、Ｋ衛星分の擬似距離誤差“εk”の共分散でなる擬似距離誤差共分散行列の逆行列として計算される。具体的には、次式（１２）で与えられる。
【数１２】
【００８３】
“μδX”は、式（８）で登場した状態変化量期待値ベクトルであり、次式（１３）で与えられる。
【数１３】
但し、“μδx，μδy，μδz，μδb”は、それぞれ“δｘ，δｙ，δｚ，δｂ”の期待値である。
【００８４】
また、“ＳδX”は、式（８）で登場した状態変化量共分散行列であり、次式（１４）で与えられる。
【数１４】
但し、４×４の各成分は、“δｘ，δｙ，δｚ，δｂ”の４成分間の共分散である。
【００８５】
また、“ＷδX”は、状態変化量逆共分散行列であり、状態変化量共分散行列“ＳδX”の逆行列として計算される。具体的には、次式（１５）で与えられる。
【数１５】
【００８６】
式（１０）を最大化する“δＸ”を求めるためには、“ｄ（ｌｏｇＬ）／ｄ（δＸ）＝０”を解けばよい。途中の計算は省略するが、最終的に、最尤推定解“δＸ”は次式（１６）のように求められる。
【数１６】
但し、“δρ”は式（４）の擬似距離変化量ベクトルであり、“Ｇ”は式（５）の幾何行列である。
【００８７】
各演算時刻それぞれについて、式（１６）に従って状態変化量ベクトル“δＸ”の最尤推定解を求める。そして、求めた状態変化量ベクトル“δＸ”の最尤推定解を用いて、状態ベクトル“Ｘ”を算出する。
【００８８】
また、本実施形態の位置算出演算では、上記の最尤推定解の算出と併せて、過去の位置算出演算の結果に基づいて状態変動正規分布モデルを更新する処理を行う。ここで、演算時刻“ｔ”における状態変化量ベクトル“δＸ”を“δＸ（ｔ）”と表記する。また、演算時刻“ｔ”における状態変化量期待値ベクトル“μδX”及び状態変化量共分散行列“ＳδX”を、それぞれ“μδX（ｔ）”及び“ＳδX（ｔ）”と表記する。
【００８９】
この場合、状態変化量期待値ベクトル“μδX（ｔ）”及び状態変化量共分散行列“ＳδX（ｔ）”を、それぞれ次式（１７）及び（１８）に従って更新する。
【数１７】
【数１８】
【００９０】
式（１７）は、１つ前の演算時刻に計算された状態変化量ベクトル“δＸ（ｔ）”で、今回の演算時刻における状態変化量期待値ベクトル“μδX（ｔ＋１）”を更新することを意味する。これは、過去の位置算出演算の結果に基づいて、演算結果（状態変化量ベクトル“δＸ”）の変動の基準（状態変化量期待値ベクトル“μδX（ｔ））を更新することに相当する。また、式（１８）は、１つ前の演算時刻に計算された状態変化量共分散行列“ＳδX（ｔ）”等の諸量を用いて、今回の演算時刻における状態変化量共分散行列“ＳδX（ｔ＋１）”を更新することを意味する。
【００９１】
状態変化量期待値ベクトル“μδX”は、状態変動正規分布モデルの期待値に係るベクトルである。また、状態変化量共分散行列“ＳδX”は、状態変動正規分布モデルの標準偏差（共分散）に係る行列である。よって、式（１７）及び（１８）は、過去の位置算出演算の結果に基づいて状態変動正規分布モデルを更新する式であると言える。
【００９２】
（２）第２の位置算出方法
第２の位置算出方法は、擬似距離誤差“ε”に正規分布モデルを適用し、状態変化量ベクトル“δＸ”の一部の算出パラメーターを確率変数として位置算出する方法である。この方法では、状態変化量ベクトル“δＸ”を構成する４つの算出パラメーター“δｘ，δｙ，δｚ，δｂ”のうちの一部の算出パラメーターを確率変数とし、残余の算出パラメーターを非確率変数として取り扱う。
【００９３】
確率変数及び非確率変数とする算出パラメーターは、適用するシステムに応じて適宜選択可能である。ここでは、ユーザーの位置変化量ベクトル“δＰ”を構成する“（δｘ，δｙ，δｚ）”の３つの算出パラメーターを確率変数とし、クロックバイアス変化量“δｂ”を非確率変数とする場合を例に挙げて説明する。この場合、式（８）の多次元状態変動正規分布モデル“ｐδX(m)（δＸ）”は、３次元状態変動正規分布モデル“ｐδX(3)（δＸ）”となる。
【００９４】
この場合、クロックバイアス変化量“δｂ”については正規分布モデルを仮定しないため、式（１０）に相当する対数尤度“ｌｏｇＬ”は、次式（１９）のように変更される。
【数１９】
【００９５】
また、式（１３）で定義した状態変化量期待値ベクトル“μδX”は、次式（２０）のように変更される。
【数２０】
但し、クロックバイアス変化量“δｂ”を非確率変数とするため、それに対応する４番目の成分をゼロとしている。
【００９６】
また、式（１４）で定義した状態変化量共分散行列“ＳδX”は、次式（２１）のように変更される。
【数２１】
但し、クロックバイアス変化量“δｂ”を非確率変数とするため、４×４の行列の４行目及び４列目の成分を全てゼロとしている。
【００９７】
この場合も、最尤推定解“δＸ”を式（１６）のように求めることができる。また、式（２０）で定義した状態変化量期待値ベクトル“μδX”と、式（２１）で定義した状態変化量共分散行列“ＳδX”とを、それぞれ式（１７）及び式（１８）に従って更新する。この場合、クロックバイアス変化量“δｂ”に対応する成分はゼロとして更新しない。
【００９８】
また、非確率変数とするクロックバイアス変化量“δｂ”については、受信信号を用いた所定の収束演算を行って解を求める。具体的には、式（３）の観測方程式に対して擬似距離誤差正規分布モデルを適用して、例えばニュートン法を利用した収束演算を行って、クロックバイアス変化量“δｂ”の近似解を求めることになる。収束演算によって算出されたクロックバイアス変化量“δｂ”は、次回の演算時刻での位置算出演算において擬似距離を近似的に求めるために使用する。
【００９９】
１−２．擬似距離誤差正規混合分布モデル（擬似距離誤差ＧＭＭ）を適用する場合
（３）第３の位置算出方法
第３の位置算出方法は、擬似距離誤差“ε”に正規混合分布モデルを適用し、状態変化量ベクトル“δＸ”の全ての算出パラメーターを確率変数として位置算出する方法である。考え方は第１の位置算出方法と同じであるが、擬似距離誤差の確率分布モデルとして、正規分布モデルではなく、正規混合分布モデルを適用する点が異なる。
【０１００】
この場合、擬似距離誤差ＧＭＭ“ｆ（ε）”と、４次元状態変動正規分布モデル“ｐδX(4)（δＸ）”とを用いて、対数尤度“ｌｏｇＬ”を次式（２２）のように定式化する。
【数２２】
【０１０１】
式（２２）において、“ｆk（εk）”は、ｋ番目の捕捉衛星の擬似距離誤差ＧＭＭを示し、次式（２３）のように一般化される。
【数２３】
【０１０２】
式（２３）において、“ｌ”は、混合（合成）する正規分布モデルの番号であり、Ｌ個（ｌ＝１，２，・・・，Ｌ）の正規分布モデルを混合する場合の式として記載している。“αl”は、ｌ番目の正規分布モデル“Ｎl”の重み（合成割合）である。“μk(l)”は、ｋ番目の捕捉衛星に係るｌ番目の正規分布モデルの期待値であり、“σk(l)”は、ｋ番目の捕捉衛星に係るｌ番目の正規分布モデルの標準偏差である。
【０１０３】
なお、上記の擬似距離誤差ＧＭＭは“ｆ（ε）”は、全ての捕捉衛星についてパラメーター値を共通化した共通のモデルとしてもよいし、捕捉衛星毎にパラメーター値を個別化した異なるモデルとしてもよい。
【０１０４】
式（２２）の対数尤度“ｌｏｇＬ”を最大化する“δＸ”を求めるためには、“ｄ（ｌｏｇＬ）／ｄ（δＸ）＝０”を解けばよい。しかし、正規混合分布（ＧＭＭ）では、対数の加算など複雑な項が存在するため、解析的に解を求めることが困難である。そこで、本実施形態では、最適化問題に対するアルゴリズムの一種であるＥＭ（Expectation Maximum）アルゴリズムを用いて最尤推定解を求める。
【０１０５】
ＥＭアルゴリズムは、確率分布モデルのパラメーターを最尤推定法に基づいて推定する手法の一種である。観測可能なデータ（サンプル）の他に、観測不可能なデータ（隠しパラメーター）が存在する場合であって、当該隠しパラメーターが確率分布モデルに依存する場合に用いられる手法である。ＥＭアルゴリズムは反復法の一種であり、尤度関数の条件付き確率に関する期待値を計算するＥ（Expectation）ステップと、期待値を最大化する解を求めるＭ（Maximum）ステップとの２つのステップを反復して最尤推定解を求める。
【０１０６】
途中の計算は省略するが、ＥＭアルゴリズムを適用すると、状態変化量ベクトル“δＸ”の最尤推定解は、次式（２４）のように求まる。
【数２４】
但し、“δＸ0”は各演算時刻における状態変化量ベクトル“δＸ”の初期値である。
【０１０７】
式（２４）における“Ｍl（δＸ0）”は、次式（２５）で与えられる。
【数２５】
【０１０８】
また、式（２５）における“ｈk（ｌ）（δＸ0）”は、次式（２６）で与えられる。
【数２６】
【０１０９】
また、式（２４）における“μl”は、次式（２７）で与えられる。
【数２７】
【０１１０】
各演算時刻において、式（２４）〜式（２７）に従って最尤推定解“δＸ”を求める。また、擬似距離誤差正規分布モデルを適用した場合と同様に、式（１３）で定義した状態変化量期待値ベクトル“μδX”と、式（１４）で定義した状態変化量共分散行列“ＳδX”とを、演算時刻毎に更新する。
【０１１１】
具体的には、状態変化量期待値ベクトル“μδX（ｔ）”及び状態変化量共分散行列“ＳδX（ｔ）”を、それぞれ次式（２８）及び（２９）に従って更新する。
【数２８】
【数２９】
【０１１２】
式（２９）における“Ｈ（δＸ）”及び“Ｆ（δＸ）”は、それぞれ次式（３０）及び（３１）で与えられる。
【数３０】
【数３１】
【０１１３】
また、式（３１）における“Ｒl”は、次式（３２）で与えられる。
【数３２】
【０１１４】
（４）第４の位置算出方法
第４の位置算出方法は、擬似距離誤差“ε”に正規混合分布モデルを適用し、状態変化量ベクトル“δＸ”の一部の算出パラメーターを確率変数として位置算出する方法である。この方法では、第２の位置算出方法と同様に、状態変化量ベクトル“δＸ”を構成する４つの算出パラメーター“δｘ，δｙ，δｚ，δｂ”のうちの一部の算出パラメーターを確率変数とし、残余の算出パラメーターを非確率変数として取り扱う。擬似距離誤差の確率分布モデルとして、正規分布モデルではなく、正規混合分布モデルを適用する点が第２の位置算出方法と異なる。
【０１１５】
この場合も、第３の位置算出方法と同様の数式を用いて、状態変化量ベクトル“δＸ”の最尤推定解を求めることができる。第３の位置算出方法と異なるのは、確率変数として取り扱う算出パラメーターについてのみ正規分布モデルを適用して計算する点である。つまり、第２の位置算出方法と同様に、状態変化量ベクトル“δＸ”の成分のうち、非確率変数として取り扱う算出パラメーターに対応する成分をゼロとして計算する。
【０１１６】
例えば、クロックバイアス変化量“δｂ”を非確率変数とする場合は、クロックバイアス変化量“δｂ”に対応する成分をゼロとして計算する。そして、クロックバイアス変化量“δｂ”については、受信信号を用いた所定の収束演算を行って解を算出する。具体的には、式（３）の観測方程式に対して擬似距離誤差正規混合分布モデル（擬似距離誤差ＧＭＭ）を適用して、例えばＥＭアルゴリズムを利用した収束演算を行って、クロックバイアス変化量“δｂ”の近似解を求めることになる。この場合、収束演算によって算出されたクロックバイアス変化量“δｂ”は、次回の演算時刻での位置算出演算において擬似距離を近似的に求めるために使用する。
【０１１７】
２．実験結果
本実施形態の位置算出方法で位置算出演算を行った実験結果について説明する。従来の手法と本実施形態の手法とのそれぞれを用いて、位置算出演算を行って実際に位置を算出する実験を行った。ユーザーに位置算出装置を携行させ、予め定められた周回経路に沿って歩行させながら位置を算出した。
【０１１８】
図９は、位置算出結果の一例を示す実験結果である。グラフの横方向は東西方向を示し、縦方向は南北方向を示す。単位はそれぞれメートル［ｍ］である。ここで示す結果は、状態変動の確率分布モデルとして正規分布を適用し、擬似距離誤差の確率分布モデルとして正規混合分布（ＧＭＭ）を適用した結果である。従来の位置算出方法で得られた軌跡を白丸のプロットで示し、本実施形態の位置算出方法で得られた軌跡を黒丸のプロットで示す。
【０１１９】
この結果を見ると、従来の位置算出方法を適用した場合は、所々に位置飛びが発生しており、それほど正確な位置は得られていないことがわかる。しかし、本実施形態の位置算出方法を適用した場合は、周回経路に沿った滑らかな軌跡が得られていることがわかる。位置算出の演算結果を独立とするのではなく、演算結果に時間的な相関を持たせて位置算出を行ったことで、算出位置の連続性が飛躍的に向上している。
【０１２０】
図１０は、上記の歩行試験において位置演算精度を計測したグラフである。グラフの横方向は東西方向の位置誤差を示し、縦方向は南北方向の位置誤差を示す。単位はそれぞれメートル［ｍ］である。グラフの中心が位置誤差“０［メートル］”に相当し、プロットがグラフの中心に近いほど位置算出精度が高いことを意味する。
【０１２１】
このグラフを見ると、従来の位置算出方法を適用した場合は、プロットが全体的に広くばらついており、位置算出精度が低くなっていることがわかる。それに対し、本実施形態の位置算出方法を適用した場合は、プロットが全体的に中心部に集中しており、従来の手法と比べて位置算出精度が格段に向上していることがわかる。
【０１２２】
図１１は、位置算出精度を示す他の実験結果の一例である。ここで示す結果は、状態変動の確率分布モデル及び擬似距離誤差の確率分布モデルとして、共に正規分布モデルを適用した結果である。図の見方は図１０と同じである。この結果を見ると、擬似距離誤差の確率分布モデルとして正規分布モデルを適用した場合も、位置誤差を示すプロットがグラフの中心部に集中しており、位置算出精度が格段に向上していることがわかる。
【０１２３】
上記の実験結果により、本実施形態の位置算出方法を用いることで、位置算出の正確性が向上することが実証された。
【０１２４】
３．実施例
次に、上記の位置算出方法を用いて位置算出を行う位置算出装置の実施例について説明する。ここでは、位置算出装置を具備した電子機器の一例である携帯型電話機の実施例として説明する。なお、本発明を適用可能な実施例が以下説明する実施例に限定されるわけでないことは勿論である。
【０１２５】
３−１．携帯型電話機の構成
図１２は、携帯型電話機１の機能構成の一例を示すブロック図である。携帯型電話機１は、ＧＰＳアンテナ５と、ＧＰＳ受信部１０と、ホスト処理部３０と、操作部４０と、表示部５０と、携帯電話用アンテナ６０と、携帯電話用無線通信回路部７０と、記憶部８０と、時計部９０とを備えて構成される。
【０１２６】
ＧＰＳアンテナ５は、ＧＰＳ衛星から発信されているＧＰＳ衛星信号を含むＲＦ（Radio Frequency）信号を受信するアンテナであり、受信信号をＧＰＳ受信部１０に出力する。
【０１２７】
ＧＰＳ受信部１０は、ＧＰＳアンテナ５から出力された信号に基づいて携帯型電話機１の位置を算出する回路或いは装置であり、ＧＰＳ受信機に相当する機能ブロックである。本実施例では、このＧＰＳ受信部１０が位置算出装置に相当する。
【０１２８】
ＧＰＳ受信部１０は、ＲＦ受信回路部１１と、ベースバンド処理回路部２０とを備えて構成される。ＧＰＳ受信部１０は、ＧＰＳ衛星信号を受信する受信部に相当する。なお、ＲＦ受信回路部１１と、ベースバンド処理回路部２０とは、それぞれ別のＬＳＩ（Large Scale Integration）として製造することも、１チップとして製造することも可能である。
【０１２９】
ＲＦ受信回路部１１は、ＲＦ信号の受信回路である。回路構成としては、例えば、ＧＰＳアンテナ５から出力されたＲＦ信号をＡ／Ｄ変換器でデジタル信号に変換し、デジタル信号を処理する受信回路を構成してもよい。また、ＧＰＳアンテナ５から出力されたＲＦ信号をアナログ信号のまま信号処理し、最終的にＡ／Ｄ変換することでデジタル信号をベースバンド処理回路部２０に出力する構成としてもよい。
【０１３０】
後者の場合には、例えば、次のようにＲＦ受信回路部１１を構成することができる。すなわち、所定の発振信号を分周或いは逓倍することで、ＲＦ信号乗算用の発振信号を生成する。そして、生成した発振信号を、ＧＰＳアンテナ５から出力されたＲＦ信号に乗算することで、ＲＦ信号を中間周波数の信号（以下、「ＩＦ（Intermediate Frequency）信号」と称す。）にダウンコンバートし、ＩＦ信号を増幅等した後、Ａ／Ｄ変換器でデジタル信号に変換して、ベースバンド処理回路部２０に出力する。
【０１３１】
ベースバンド処理回路部２０は、ＲＦ受信回路部１１から出力された受信信号に対して、キャリア除去や相関演算等を行ってＧＰＳ衛星信号を捕捉する。そして、捕捉したＧＰＳ衛星信号から抽出した時刻情報や衛星軌道情報等を利用して、携帯型電話機１の位置及び時計誤差を算出する。
【０１３２】
ホスト処理部３０は、記憶部８０に記憶されているシステムプログラム等の各種プログラムに従って携帯型電話機１の各部を統括的に制御するプロセッサーであり、ＣＰＵ（Central Processing Unit）等のプロセッサーを有して構成される。ホスト処理部３０は、ベースバンド処理回路部２０から取得した位置座標をもとに、表示部５０に現在位置を指し示した地図を表示させたり、その位置座標を各種のアプリケーション処理に利用する。
【０１３３】
操作部４０は、例えばタッチパネルやボタンスイッチ等により構成される入力装置であり、押下されたキーやボタンの信号をホスト処理部３０に出力する。この操作部４０の操作により、通話要求やメール送受信要求、位置算出要求等の各種指示入力がなされる。
【０１３４】
表示部５０は、ＬＣＤ（Liquid Crystal Display）等により構成され、ホスト処理部３０から入力される表示信号に基づいた各種表示を行う表示装置である。表示部５０には、位置表示画面や時刻情報等が表示される。
【０１３５】
携帯電話用アンテナ６０は、携帯型電話機１の通信サービス事業者が設置した無線基地局との間で携帯電話用無線信号の送受信を行うアンテナである。
【０１３６】
携帯電話用無線通信回路部７０は、ＲＦ変換回路、ベースバンド処理回路等によって構成される携帯電話の通信回路部であり、携帯電話用無線信号の変調・復調等を行うことで、通話やメールの送受信等を実現する。
【０１３７】
記憶部８０は、ＲＯＭ（Read Only Memory）やフラッシュＲＯＭ、ＲＡＭ（Random Access Memory）等の記憶装置を有して構成され、ホスト処理部３０が携帯型電話機１を制御するためのシステムプログラムや、各種アプリケーション処理を実行するための各種プログラムやデータ等を記憶する。
【０１３８】
時計部９０は、携帯型電話機１の内部時計であり、水晶振動子及び発振回路でなる水晶発振器等を有して構成される。時計部９０の計時時刻は、ベースバンド処理回路部２０及びホスト処理部３０に随時出力される。時計部９０は、ベースバンド処理回路部２０により算出された時計誤差を用いて較正される。
【０１３９】
３−２．ベースバンド処理回路部の構成
図１３は、ベースバンド処理回路部２０の回路構成及び記憶部２３のデータ構成の説明図である。ベースバンド処理回路部２０は、主要な機能構成として、処理部２１と、記憶部２３とを備える。
【０１４０】
処理部２１は、ベースバンド処理回路部２０の各機能部を統括的に制御する制御装置及び演算装置であり、ＣＰＵやＤＳＰ（Digital Signal Processor）等のプロセッサーを有して構成される。処理部２１は、衛星捕捉部２１１と、位置算出部２１３とを機能部として有する。
【０１４１】
衛星捕捉部２１１は、ＧＰＳ衛星の捕捉を行う機能部である。具体的には、ＲＦ受信回路部１１から出力されるデジタル化された受信信号に対して、キャリア除去及び相関演算のデジタル信号処理を実行して、ＧＰＳ衛星を捕捉する。そして、相関演算結果に対するピーク判定を行い、受信キャリア信号のドップラー周波数や受信周波数、受信Ｃ／Ａコードのコード位相といった情報をメジャメント情報として取得する。
【０１４２】
また、衛星捕捉部２１１は、相関演算結果に基づいて航法データを復号する。受信キャリア信号の位相（キャリア位相）と、受信Ｃ／Ａコードの位相（コード位相）とが検出され、相関がとれた状態になると、相関値の時間変化をもとに、航法データを構成する各ビット値をデコードすることができる。この位相同期は、例えば位相ロックループとして知られるＰＬＬ（Phase Locked Loop）処理により実現される。
【０１４３】
位置算出部２１３は、衛星捕捉部２１１により捕捉された各ＧＰＳ衛星に係るメジャメント情報や航法データ、時刻情報、衛星情報といった諸量を用いて、上記の原理に従って位置算出演算を行って、携帯型電話機１の位置（位置座標）及び時計誤差（クロックバイアス）を算出する。そして、算出した位置をホスト処理部３０に出力するとともに、算出した時計誤差で時計部９０を較正する。
【０１４４】
位置算出部２１３は、ＧＰＳ受信部１０による受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う位置算出部に相当する。
【０１４５】
記憶部２３は、ベースバンド処理回路部２０のシステムプログラムや、衛星捕捉機能、位置算出機能といった各種機能を実現するための各種プログラム、データ等を記憶する。また、各種処理の処理中データ、処理結果などを一時的に記憶するワークエリアを有する。
【０１４６】
記憶部２３には、プログラムとして、処理部２１により読み出され、ベースバンド処理（図１４参照）として実行される第１のベースバンド処理プログラム２３１が記憶されている。また、ベースバンド処理プログラム２３１は、各種の位置算出処理（図１５〜図１８参照）として実行される位置算出プログラム２３１１をサブルーチンとして含む。原理で説明した第１〜第４の位置算出方法の何れかに対応するプログラムが位置算出プログラム２３１１として記憶されている。
【０１４７】
また、記憶部２３には、主要なデータとして、衛星軌道データ２３２と、捕捉した各ＧＰＳ衛星に係る衛星別データ２３３と、状態変動正規分布モデル２３４と、演算結果データ２３５とが記憶される。
【０１４８】
衛星軌道データ２３２は、全てのＧＰＳ衛星の概略の衛星軌道情報を記憶したアルマナックや、各ＧＰＳ衛星それぞれについて詳細な衛星軌道情報を記憶したエフェメリス等のデータである。この衛星軌道データ２３２は、ＧＰＳ衛星から受信したＧＰＳ衛星信号をデコードすることで取得する他、例えば携帯型電話機１の基地局やアシストサーバーからアシストデータとして取得する。
【０１４９】
衛星別データ２３３は、メジャメント情報２３３Ａと、衛星情報２３３Ｂと、測位計算用諸量２３３Ｃとを含む。メジャメント情報２３３Ａは、当該捕捉衛星について測定したＧＰＳ衛星信号のコード位相や受信周波数、ドップラー周波数といった情報である。衛星情報２３３Ｂは、当該捕捉衛星の位置や速度、移動方向といった情報である。測位計算用諸量２３３Ｃは、上記の原理に従って測位計算を行う際に用いる各種の諸量である。
【０１５０】
状態変動正規分布モデル２３４は、図１で図示・説明した状態変動正規分布モデルであり、モデルパラメーター値として、例えば、状態変化量期待値２３４Ａと、状態変化量標準偏差２３４Ｂとを含む。状態変動正規分布モデルは、演算時刻毎に更新される。
【０１５１】
演算結果データ２３５は、位置算出演算の結果として得られたデータであり、携帯型電話機１の位置ベクトル“Ｐ”やクロックバイアス“ｂ”がこれに含まれる。
【０１５２】
３−３．処理の流れ
図１４は、記憶部２３に記憶されているベースバンド処理プログラム２３１が処理部２１により読み出されることで、ベースバンド処理回路部２０において実行されるベースバンド処理の流れを示すフローチャートである。
【０１５３】
最初に、処理部２１は、記憶部２３にサブルーチンとして記憶されている位置算出プログラム２３１１に従って位置算出処理を開始する（ステップＡ１）。位置算出処理は、図１５〜図１８の第１〜第４の位置算出処理の何れかの処理である。
【０１５４】
次いで、処理部２１は、演算結果の出力タイミングであるか否かを判定する（ステップＡ３）。演算結果の出力タイミングは、任意のタイミングを設定可能である。例えば、位置算出演算の時間間隔と同じ時間間隔のタイミングとしてもよいし、位置算出演算の時間間隔よりも長い時間間隔のタイミングとしてもよい。また、ユーザーから演算結果の出力が指示されたタイミングを出力タイミングとしてもよい。
【０１５５】
出力タイミングであると判定した場合は（ステップＡ３；Ｙｅｓ）、処理部２１は、最新の演算結果をホスト処理部３０に出力する（ステップＡ５）。ステップＡ３において出力タイミングではないと判定した場合（ステップＡ３；Ｎｏ）、又は、ステップＡ５の後、処理部２１は、ベースバンド処理を終了するか否かを判定する（ステップＡ７）。
【０１５６】
ベースバンド処理を終了しないと判定した場合は（ステップＡ７；Ｎｏ）、処理部２１は、ステップＡ３に戻る。また、終了すると判定した場合は（ステップＡ７；Ｙｅｓ）、処理部２１は、ベースバンド処理を終了する。
【０１５７】
図１５は、第１の位置算出方法に対応する第１の位置算出処理の流れを示すフローチャートである。最初に、処理部２１は、状態ベクトル“Ｘ”の初期値を設定する（ステップＢ１）。具体的には、サーバーアシストによって携帯型電話機１の基地局から取得した位置（基地局位置）や、予め定められた位置（例えば固定値）を、初期位置として設定する。また、クロックバイアス“ｂ”にも所定の初期値（例えば固定値）を設定する。
【０１５８】
次いで、処理部２１は、捕捉対象衛星選定処理を行う（ステップＢ３）。具体的には、時計部９０で計時されている現在時刻において、初期位置の天空に位置するＧＰＳ衛星を、記憶部２３に記憶されたアルマナックやエフェメリス等の衛星軌道データ２３２を用いて判定して、捕捉対象衛星に選定する。
【０１５９】
次いで、処理部２１は、状態変化量推定処理を行う（ステップＢ５〜Ｂ１７）。具体的には、処理部２１は、各捕捉対象衛星についてループＡを実行する（ステップＢ５〜Ｂ１３）。ループＡでは、処理部２１は、当該捕捉対象衛星からＧＰＳ衛星信号を受信して擬似距離を測定し、その値を擬似距離測定値“ρc”とする（ステップＢ７）。擬似距離測定値“ρc”は、位相サーチの結果として得られるコード位相を用いて算出することができる。
【０１６０】
また、処理部２１は、擬似距離近似値“ρa”を算出する（ステップＢ９）。具体的には、衛星軌道データ２３２から求まる当該捕捉対象衛星の位置ベクトル“Ｐk”と、携帯型電話機１の最新の位置ベクトル“Ｐ”及び最新のクロックバイアス“ｂ”とを用いて、擬似距離近似値を“ρa＝‖Ｐk−Ｐ‖＋ｂ”として算出する。但し、この式では、クロックバイアス“ｂ”の単位を距離に換算している。
【０１６１】
その後、処理部２１は、擬似距離変化量“δρ”を算出する（ステップＢ１１）。具体的には、ステップＢ７で測定した擬似距離測定値“ρc”と、ステップＢ９で算出した擬似距離近似値“ρa”との差を算出して、擬似距離変化量“δρ”とする。そして、処理部２１は、次の捕捉対象衛星へと処理を移行する。全ての捕捉対象衛星についてステップＢ７〜Ｂ１１を行った後、処理部２１は、ループＡを終了する（ステップＢ１３）。
【０１６２】
次いで、処理部２１は、幾何行列Ｇを算出する（ステップＢ１５）。具体的には、衛星軌道データ２３２から求まる各捕捉対象衛星の衛星位置ベクトルと、携帯型電話機１の最新の位置ベクトルとを用いて、携帯型電話機１から各捕捉対象衛星に向かう視線ベクトル“ｌ”を計算する。そして、式（５）の幾何行列Ｇを求める。
【０１６３】
その後、処理部２１は、擬似距離誤差正規分布モデル及び状態変動正規分布モデルを適用して、状態変化量“δＸ”の最尤推定解を求める（ステップＢ１７）。具体的には、式（１６）に従って状態変化量“δＸ”を計算する。擬似距離誤差“ε”に正規分布を仮定するため、ここで得られる最尤推定解は最小二乗解に相当する。そして、処理部２１は、状態変化量推定処理を終了する。
【０１６４】
次いで、処理部２１は、ステップＢ１７で最尤推定解として推定した状態変化量ベクトル“δＸ”を用いて状態ベクトル“Ｘ”を修正する（ステップＢ１９）。具体的には、最尤推定解として求めた状態変化量ベクトル“δＸ”を最新の状態ベクトル“Ｘ”に加算することで、今回の演算時刻における状態ベクトル“Ｘ”を新たに算出する。
【０１６５】
その後、処理部２１は、位置算出を終了するか否かを判定し（ステップＢ２３）、まだ終了しないと判定した場合は（ステップＢ２３；Ｎｏ）、状態変動正規分布モデル更新処理を行う（ステップＢ２５）。具体的には、状態変化量期待値ベクトル“μδX”及び状態変化量共分散行列“ＳδX”を、それぞれ式（１７）及び式（１８）に従って更新する。
【０１６６】
状態変動正規分布モデル更新処理を行った後、処理部２１は、ステップＢ３に戻り、次の演算時刻での位置算出演算へと移行する。また、ステップＢ２３において位置算出を終了すると判定した場合は（ステップＢ２３；Ｙｅｓ）、処理部２１は、第１の位置算出処理を終了する。
【０１６７】
ステップＢ３〜Ｂ２５までの一連の処理が、一の演算時刻での位置算出演算に相当する。第１の位置算出処理において特徴的であるのは、各演算時刻での位置算出演算において収束演算が省略されている点である。従来の位置算出方法では、最小二乗法を利用して解を求める際に、例えばニュートン法を利用した収束演算を行って近似解を求めていた。しかし、本実施形態では、演算結果の変動を確率変数とし、その上、演算時刻毎に確率分布モデルを更新することにしているため、収束演算を不要とする。
【０１６８】
図１６は、第２の位置算出方法に対応する第２の位置算出処理の流れを示すフローチャートである。状態変化量ベクトル“δＸ”のうち、位置変化量ベクトル“δＰ”の成分“δｘ，δｙ，δｚ”を確率変数とし、クロックバイアス変化量“δｂ”を非確率変数として位置算出を行うフローチャートである。なお、図１５の第１の位置算出処理と同一のステップについては同一の符号を付して説明を省略する。
【０１６９】
処理部２１は、ステップＢ３において捕捉対象衛星選定処理を行った後、状態変化量推定処理を行う（ステップＣ５）。具体的には、処理部２１は、第１の位置算出処理の状態変化量推定処理（ステップＢ５〜Ｂ１７）と同様の処理を行う。但し、ここではクロックバイアス“δｂ”を非確率変数とするため、状態変化量ベクトル“δＸ”のクロックバイアス変化量“δｂ”に対応する成分をゼロとして計算する。つまり、状態変化量推定処理は、位置変化量ベクトル“δＰ”を推定する位置変化量推定処理となる。
【０１７０】
また、状態変化量推定処理では、１回前の位置算出演算で収束演算によって算出したクロックバイアス“ｂ”の値を用いて、擬似距離近似値“ρa”を算出する（ステップＢ９）。具体的には、衛星軌道データ２３２から求まる当該捕捉対象衛星の位置ベクトル“Ｐk”と、前回の位置算出演算で算出した位置ベクトル“Ｐ”と、前回の位置算出演算で収束演算によって算出したクロックバイアス“ｂ”とを用いて、擬似距離近似値“ρa＝‖Ｐk−Ｐ‖＋ｂ”を算出する。これは、過去に収束演算によって算出された所定のパラメーターの値を用いて位置算出演算を行うことに相当する。
【０１７１】
次いで、処理部２１は、クロックバイアス変化量推定処理を行う（ステップＣ７〜Ｃ１９）。クロックバイアス変化量推定処理は、ＧＰＳ衛星信号の受信信号を用いた所定の収束演算を行って、非確率変数とするクロックバイアス変化量“δｂ”の近似解を求める処理である。
【０１７２】
詳細には、処理部２１は、各捕捉対象衛星についてループＢを行う（ステップＣ７〜Ｃ１３）。ループＢでは、処理部２１は、擬似距離近似値“ρa”を算出する（ステップＣ９）。そして、処理部２１は、擬似距離変化量“δρ”を算出する（ステップＣ１１）。具体的には、状態変化量推定処理のステップＢ７で測定した擬似距離測定値“ρc”と、ステップＣ９で算出した擬似距離近似値“ρa”との差を算出して、擬似距離変化量“δρ”とする。
【０１７３】
その後、処理部２１は、次の捕捉対象衛星へと処理を移行する。全ての捕捉対象衛星についてステップＣ９及びＣ１１を行った後、処理部２１は、ループＢを終了する（ステップＣ１３）。
【０１７４】
次いで、処理部２１は、擬似距離誤差正規分布モデルを適用して、クロックバイアス変化量“δｂ”の最小二乗解を求める（ステップＣ１７）。そして、処理部２１は、求めたクロックバイアス変化量“δｂ”について、所定の収束条件が成立したか否かを判定する（ステップＣ１８）。具体的には、例えば、今回の反復ステップで求めたクロックバイアス変化量“δｂ”と、１回前の反復ステップで求めたクロックバイアス変化量“δｂ”との差の絶対値が所定の閾値よりも小さくなった場合に、収束条件が成立したと判定する。
【０１７５】
収束条件が成立していないと判定した場合は（ステップＣ１８；Ｎｏ）、処理部２１は、ステップＣ１７で求めたクロックバイアス変化量“δｂ”を用いて、クロックバイアス“ｂ”を修正して更新する（ステップＣ１９）。そして、処理部２１は、ステップＣ７に戻る。また、収束条件が成立したと判定した場合は（ステップＣ１８；Ｙｅｓ）、処理部２１は、クロックバイアス変化量推定処理を終了する。
【０１７６】
その後、処理部２１は、状態変化量推定処理において推定した位置変化量ベクトル“δＰ”と、クロックバイアス変化量推定処理において推定したクロックバイアス変化量“δｂ”とを状態変化量ベクトル“δＸ”とし、これを用いて状態ベクトル“Ｘ”を修正する（ステップＣ２１）。つまり、最新の状態ベクトル“Ｘ”に状態変化量ベクトル“δＸ”を加算して、今回の演算時刻における状態ベクトル“Ｘ”を新たに算出する。
【０１７７】
その後、処理部２１は、位置算出を終了するか否かを判定し（ステップＢ２３）、まだ終了しないと判定した場合は（ステップＢ２３；Ｎｏ）、状態変動正規分布モデル更新処理を行う（ステップＢ２５）。その後、処理部２１は、ステップＢ３に戻り、次の演算時刻での位置算出演算へと移行する。また、ステップＢ２３において位置算出を終了すると判定した場合は（ステップＢ２３；Ｙｅｓ）、処理部２１は、第２の位置算出処理を終了する。
【０１７８】
この場合も、ステップＢ３〜Ｂ２５までの一連の処理が、一の演算時刻での位置算出演算に相当する。また、ステップＣ７〜Ｃ１９の繰り返し処理が、受信信号を用いた所定の収束演算に相当する。第２の位置算出処理において特徴的であるのは、確率変数として取り扱う位置変化量ベクトル“δＰ”については収束演算を省略しているが、非確率変数として取り扱うクロックバイアス変化量“δｂ”については収束演算を省略していない点である。クロックバイアス変化量“δｂ”は確率変数とせず、全くの未知数として取り扱うため、収束演算を行わなければ最適解が求まらないためである。
【０１７９】
図１７は、第３の位置算出方法に対応する第３の位置算出処理の流れを示すフローチャートである。第３の位置算出処理の基本的な流れは、図１５の第１の位置算出処理と同じである。しかし、第１の位置算出処理のステップＢ１７がステップＤ１７に置き換わる点が大きく異なる。
【０１８０】
ステップＤ１７では、処理部２１は、擬似距離誤差正規混合分布モデル及び状態変動正規分布モデルを適用して、状態変化量ベクトル“δＸ”の最尤推定解を求める。具体的には、処理部２１は、式（２４）〜式（２７）に従って、状態変化量ベクトル“δＸ”の最尤推定解を計算する。このステップは、ＥＭアルゴリズムを適用した計算ステップである。
【０１８１】
また、第３の位置算出処理では、ステップＢ２５の状態変動正規分布モデル更新処理において、処理部２１は、状態変化量期待値ベクトル“μδX”を式（２８）に従って更新する。また、処理部２１は、状態変化量共分散行列“ＳδX”を式（２９）〜式（３２）に従って更新する。
【０１８２】
図１８は、第４の位置算出方法に対応する第４の位置算出処理の流れを示すフローチャートである。第４の位置算出処理の基本的な流れは、図１６の第２の位置算出処理と同じである。しかし、第２の位置算出処理のステップＣ５がステップＥ５に置き換わり、ステップＣ１７がステップＥ１７に置き換わる点が大きく異なる。
【０１８３】
第４の位置算出処理では、ステップＥ５の状態変化量推定処理において、クロックバイアス“δｂ”を非確率変数として取り扱い、状態変化量ベクトル“δＸ”のクロックバイアス変化量“δｂ”に対応する成分をゼロとして計算を行う。つまり、この場合の状態変化量推定処理は、位置変化量ベクトル“δＰ”を推定する位置変化量推定処理となる。
【０１８４】
また、第２の位置算出処理と同様に、状態変化量推定処理では、１回前の位置算出演算で収束演算によって算出されたクロックバイアス“ｂ”の値を用いて、擬似距離近似値“ρa”を算出する（ステップＢ９）。つまり、過去に収束演算によって算出された所定のパラメーターの値を用いて位置算出演算を行う。
【０１８５】
また、第４の位置算出処理では、ステップＣ７〜ステップＣ１９のクロックバイアス変化量推定処理において、擬似距離誤差正規混合分布モデルを適用して、例えばＥＭアルゴリズムを利用した収束演算を行って（ステップＣ１８；Ｎｏ→ステップＣ１９→ステップＣ７）、クロックバイアス変化量“δｂ”の最尤推定解を求める（ステップＥ１７）。
【０１８６】
そして、処理部２１は、クロックバイアス変化量“δｂ”について所定の収束条件が成立したか否かを判定し（ステップＣ１８）、成立していないと判定した場合は（ステップＣ１８；Ｎｏ）、ステップＥ１７で求めたクロックバイアス変化量“δｂ”を用いてクロックバイアス“ｂ”を修正・更新した後（ステップＣ１９）、ステップＣ７に戻る。また、収束条件が成立したと判定した場合は（ステップＣ１８；Ｙｅｓ）、処理部２１は、ステップＣ２１へと移行する。
【０１８７】
４．作用効果
本実施形態によれば、測位用衛星信号の一種であるＧＰＳ衛星信号を受信する。そして、ＧＰＳ衛星信号の受信信号を用いた演算であって、ＧＰＳ受信機の位置及びクロックバイアスの変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う。
【０１８８】
具体的には、ＧＰＳ受信機（位置算出装置）の位置及びクロックバイアスの変動を確率変数とする状態変動正規分布モデルを設定する。また、擬似距離に含まれる誤差の変動を確率変数とする確率分布モデルとして、擬似距離誤差正規分布モデル又は擬似距離誤差正規混合分布モデルを設定する。そして、これらの確率モデルを用いて、最尤推定法に基づく演算を行って、位置及びクロックバイアスの変化量の最尤推定解を求める。これにより、演算時刻毎の演算結果を独立とするのではなく、演算結果に時間的な相関を持たせて位置算出を行うことが可能となる。この場合、測定した擬似距離に瞬間的に大きな誤差が混入したとしても、位置算出精度が低下することを防止できる。
【０１８９】
また、本実施形態では、過去の位置算出演算の結果に基づいて演算結果の変動の基準を更新し、更新した基準を用いて演算結果の変動を算出する。詳細には、１つ前の演算時刻に算出された位置及びクロックバイアスの変化量を用いて、今回の演算時刻での位置及びクロックバイアスの変化量の基準を更新する。これは、状態変動正規分布モデルの期待値を更新することに相当する。
【０１９０】
また、本実施形態では、１つ前の演算時刻での位置及びクロックバイアスの変化量の共分散を用いて、今回の演算時刻での位置及びクロックバイアスの変化量の共分散を更新する。これは、状態変動正規分布モデルの標準偏差を更新することに相当する。つまり、演算結果の変動に係る正規分布モデルを固定的に設定するのではなく、１つ前の演算時刻の演算結果に基づいて当該正規分布モデルを更新する。これにより、位置及びクロックバイアスの時間変化に俊敏に対応することができるようになり、移動体の位置をより正確に求めることが可能となる。
【０１９１】
５．変形例
５−１．ドップラー測位
上記の実施形態では、擬似距離を用いた位置算出方法を例に挙げて説明したが、ドップラーを利用した位置算出方法に対しても、本発明を実質的に同一に適用可能である。ドップラーを利用した位置算出方法は「ドップラー測位」と呼ばれる。
【０１９２】
ドップラーに関して、次式（３３）のモデル式が成立する。
【数３３】
但し、単位を全て速度の次元［ｍ／ｓ］として式を記載している。
【０１９３】
式（３３）において、“ｈk”は、ｋ番目の衛星について測定した周波数のずれ（以下、「周波数ずれ」と称す。）である。つまり、ＧＰＳ衛星信号の受信信号の周波数（受信周波数）と規定された搬送波周波数（１．５７５４２［ＧＨｚ］）との差として計算される。受信周波数は、周波数方向の相関演算を行うことで取得することができる。
【０１９４】
“Ｖk”は、ｋ番目の捕捉衛星の速度ベクトルであり、“Ｖ”は、ＧＰＳ受信機の速度ベクトルである。従って、“ＶR＝Ｖk−Ｖ”は、ＧＰＳ受信機とｋ番目の捕捉衛星との相対速度ベクトルである。また、“ｌk”は、ＧＰＳ受信機からｋ番目の捕捉衛星に向かう視線ベクトルである。相対速度ベクトル“ＶR＝Ｖk−Ｖ”を視線ベクトル“ｌk”の方向に投影することで、視線方向の相対速度ベクトル（ドップラーシフト）が求まる。
【０１９５】
また、式（３３）において、“ｄ”は、ＧＰＳ受信機の動作クロックのドリフト（クロックドリフト）である。また、“εkh”は、ｋ番目の衛星の周波数ずれに内在する観測誤差である。
【０１９６】
式（３３）に対してニュートン法を適用することを考える。ＧＰＳ受信機の位置ベクトル“Ｐ”、クロックバイアス“ｂ”、速度ベクトル“Ｖ”及びクロックドリフト“ｄ”を成分とする状態ベクトル“Ｘ＝（Ｐ，ｂ，Ｖ，ｄ）T”を定義する。つまり、位置情報の他に、クロックバイアス、速度情報及びクロックドリフトを状態ベクトル“Ｘ”の成分とする。
【０１９７】
ここで、式（３３）の（右辺）−（左辺）で表される次式（３４）の関数“ｇ（Ｐ，Ｖ）”を定義する。
【数３４】
【０１９８】
関数“ｇ”は、速度ベクトル“Ｖ”及び視線ベクトル“ｌ”を含む。視線ベクトル“ｌ”は、位置ベクトル“Ｐ”を変数とする。そのため、関数“ｇ”は、位置ベクトル“Ｐ”及び速度ベクトル“Ｖ”を変数とする関数である。
【０１９９】
関数“ｇ（Ｐ，Ｖ）”の位置ベクトル“Ｐ”、速度ベクトル“Ｖ”、クロックバイアス“ｂ”及びクロックドリフト“ｄ”での偏微分を求めると、それぞれ次式（３５）〜（３８）のようになる。
【数３５】
【数３６】
【数３７】
【数３８】
【０２００】
ここで、捕捉したＫ個の衛星について関数“ｇ（Ｐ，Ｖ）”を纏めた関数“Ｇ（Ｐ，Ｖ）”を便宜的に観測関数と呼び、次式（３９）のように定義する。
【数３９】
【０２０１】
この場合、次式（４０）の観測方程式が成立する。
【数４０】
【０２０２】
式（４０）の左辺は、位置ベクトルの初期値“Ｐ0”と、クロックバイアスの初期値“ｂ0”と、速度ベクトルの初期値“Ｖ0”と、クロックドリフトの初期値“ｄ0”とにより定まる観測関数の値“Ｇ（Ｐ0，ｂ0，Ｖ0，ｄ0）”である。
【０２０３】
また、式（４０）の右辺において、ベクトル“（δＰ，δｂ，δＶ，δｄ）T”を状態ベクトル“Ｘ”の変動を示す状態変化量ベクトル“δＸ”と定義する。また、便宜的に、観測関数“Ｇ”の各変数での偏微分で定義される行列を偏微分行列“Ｙ”と定義する。状態変化量ベクトル“δＸ”には、各変数の初期値により定まる偏微分行列“Ｙ”の値が乗算される。
【０２０４】
このように、式（３３）から、状態変化量ベクトル“δＸ＝（δＰ，δｂ，δＶ，δｄ）T”を未知数とする式（４０）の観測方程式が導出される。この観測方程式を利用して、上記の実施形態と同様に、状態変化量ベクトル“δＸ”の最尤推定解を求める。上記の実施形態と異なるのは、位置算出演算の基礎となる観測方程式が、式（３）から式（４０）に置き換わる点である。
【０２０５】
この場合、状態変化量ベクトル“δＸ＝（δＰ，δｂ，δＶ，δｄ）T”に含まれる算出パラメーター“（δｘ，δｙ，δｚ，δｂ，δｖx，δｖy，δｖz，δｄ）”の全部又は一部を確率変数とする確率分布モデルを定義する。位置情報の他に、クロックバイアス、速度情報及びクロックドリフトの変動を確率変数として取り扱う。
【０２０６】
８つの算出パラメーターの全てを確率変数とする場合は、式（８）の多次元状態変動正規分布モデルを、８次元状態変動正規分布モデル“ｐδX(8)（δＸ）”として定義する。また、８つの算出パラメーターのうち、例えばクロックバイアス“ｂ”及びクロックドリフト“ｄ”を除いた６つの算出パラメーターを確率変数とするのであれば、式（８）の多次元状態変動正規分布モデルを、６次元状態変動正規分布モデル“ｐδX(6)（δＸ）”として定義する。
【０２０７】
なお、式（４０）では未知数が８個であるため、解を求めるためには、少なくとも８個の衛星を捕捉する必要がある。そこで、未知数を減らす工夫をすると実用的である。例えば、速度ベクトル“Ｖ”を固定値（例えばゼロ）とし、位置ベクトル“Ｐ”、クロックバイアス“ｂ”及びクロックドリフト“ｄ”を未知数とする。この場合、式（４０）の観測方程式は、次式（４１）のように書き換えることができる。
【数４１】
【０２０８】
式（４１）では、状態ベクトル“Ｘ”の変動を示す状態変化量ベクトルを“δＸ＝（δＰ，δｂ，δｄ）T”と定義している。速度ベクトル“Ｖ＝（ｖx，ｖy，ｖz）”を固定値とすることで、未知数を８個から５個まで減らすことができる。これにより、少なくとも５個の衛星を捕捉することができれば、式（４１）の解を求めることができる。
【０２０９】
図１９は、上記の実施例における処理部２１が実行するドップラー測位処理の流れを示すフローチャートである。このドップラー測位処理は、図１４のベースバンド処理のステップＡ１において位置算出処理に置き換えて適用可能な処理である。また、このドップラー測位処理は、式（４０）の観測方程式に基づく処理であり、状態変化量ベクトル“δＸ”の全ての算出パラメーターを確率変数とする図１５の第１の位置算出処理に対応する処理である。
【０２１０】
最初に、処理部２１は、状態ベクトル“Ｘ”の初期値を設定する（ステップＦ１）。つまり、位置ベクトル“Ｐ”、クロックバイアス“ｂ”、速度ベクトル“Ｖ”及びクロックドリフト“ｄ”に所定の初期値を設定する。
【０２１１】
次いで、処理部２１は、捕捉対象衛星選定処理を行う（ステップＦ３）。そして、処理部２１は、状態変化量推定処理を行う（ステップＦ５〜Ｆ１９）。具体的には、処理部２１は、各捕捉対象衛星についてループＣを実行する（ステップＦ５〜Ｆ１３）。
【０２１２】
ループＣでは、処理部２１は、当該捕捉対象衛星からＧＰＳ衛星信号を受信して周波数ずれ“ｈ”を算出する（ステップＦ７）。周波数ずれ“ｈ”は、周波数サーチを行うことで取得した受信信号の周波数（受信周波数）と、規定された搬送波の周波数（規定搬送波周波数）との差として算出される。この周波数ずれ“ｈ”は、ドップラー測位における観測量となる。
【０２１３】
次いで、処理部２１は、位置算出装置と当該捕捉対象衛星との相対速度ベクトル“ＶR”を算出する（ステップＦ９）。相対速度ベクトル“ＶR”は、位置算出装置の最新の速度ベクトル“Ｖ”と、当該捕捉対象衛星の最新の速度ベクトル“Ｖk”とを用いて算出される。初回測位の場合は、速度ベクトルの初期値“Ｖ0”を最新の速度ベクトルとし、２回目以降の測位の場合は、前回の測位で算出した速度ベクトルを最新の速度ベクトルとして計算する。
【０２１４】
その後、処理部２１は、位置算出装置から当該捕捉対象衛星に向かう視線ベクトル“ｌ”を算出する（ステップＦ１１）。視線ベクトル“ｌ”は、位置算出装置の最新の位置ベクトル“Ｐ”と、当該捕捉対象衛星の最新の位置ベクトル“Ｐk”とを用いて算出される。初回測位の場合は、位置ベクトルの初期値“Ｐ0”を最新の位置ベクトルとし、２回目以降の測位の場合は、前回の測位で算出した位置ベクトルを最新の位置ベクトルとして計算する。
【０２１５】
その後、処理部２１は、次の捕捉対象衛星へと処理を移行する。全ての捕捉対象衛星についてステップＦ７〜Ｆ１１を行った後、処理部２１は、ループＣを終了する（ステップＦ１３）。
【０２１６】
ループＣの後、処理部２１は、最新の位置ベクトル“Ｐ”、クロックバイアス“ｂ”、速度ベクトル“Ｖ”及びクロックドリフト“ｄ”を用いて、式（４０）の左辺の観測関数“Ｇ”の値を算出する（ステップＦ１５）。また、処理部２１は、最新の位置ベクトル“Ｐ”、クロックバイアス“ｂ”、速度ベクトル“Ｖ”及びクロックドリフト“ｄ”を用いて、式（４０）の右辺の偏微分行列“Ｙ”の値を算出する（ステップＦ１７）。
【０２１７】
次いで、処理部２１は、式（４０）の観測方程式に従って、観測誤差正規分布モデル及び状態変動正規分布モデルを適用して、状態変化量“δＸ”の最尤推定解を求める（ステップＦ１９）。観測誤差正規分布モデルは、観測量である周波数ずれ“ｈ”に含まれる観測誤差“εh”を確率変数とし、その分布として正規分布を仮定したモデルである。
【０２１８】
ステップＦ１９を行った後、処理部２１は、状態変化量推定処理を終了する。以降のステップは、第１の位置算出処理と同様である。
【０２１９】
なお、ここではフローチャートの図示及び説明を省略するが、図１６〜図１８の第２〜第４の位置算出処理に対応するドップラー測位処理についても、それぞれ同様に構成することが可能である。
【０２２０】
５−２．確率変数
上記の実施形態では、位置算出演算の演算結果の変動を確率変数とする確率分布モデルを定義したが、位置算出演算の演算結果そのものを確率変数とする確率分布モデルを定義して位置算出演算を行うことも可能である。
【０２２１】
位置ベクトル“Ｐ”に着目して考えてみる。位置ベクトル“Ｐ”は、位置ベクトルの初期値“Ｐ0”に位置変化量ベクトル“δＰ”を加算することで求められる。位置ベクトルの初期値“Ｐ0”は固定値であるため、位置変化量ベクトル“δＰ＝Ｐ−Ｐ0”を確率変数として取り扱うことと、位置ベクトル“Ｐ”そのものを確率変数として取り扱うこととは同値である。従って、位置変化量ベクトル“δＰ”の代わりに位置ベクトル“Ｐ”そのものを確率変数として位置算出演算を行うことも可能である。クロックバイアス“ｂ”、速度ベクトル“Ｖ”及びクロックドリフト“ｄ”についても同様である。
【０２２２】
この場合は、例えばＧＰＳ受信機の状態ベクトル“Ｘ”そのものを確率変数とする状態正規分布モデルを定義する。そして、当該状態正規分布モデルを用いて、上記の実施形態と同様に位置算出演算を行う。状態正規分布モデルの更新では、過去に算出された位置ベクトル“Ｐ”やクロックバイアス“ｂ”を用いて、状態正規分布モデルの期待値や標準偏差を更新する。
【０２２３】
５−３．状態変動正規分布モデルの更新
上記の実施形態では、１時刻前の位置算出演算の結果に基づいて状態変動正規分布モデルを更新したが、あくまでもこれは一例である。例えば、２時刻前や３時刻前の位置算出演算の結果に基づいて状態変動正規分布モデルを更新してもよい。
【０２２４】
また、過去所定期間分（例えば５時刻分や１０時刻分）の位置算出演算の結果を平均処理し、その結果に基づいて状態変動正規分布モデルを更新してもよい。平均処理としては、単純な算術平均や幾何平均を適用してもよいし、加重平均を適用してもよい。加重平均を適用する場合は、現在の演算時刻により近い時刻での位置算出演算の結果ほど、重みをより高く設定して加重平均すると効果的である。
【０２２５】
また、外部から取得可能な参照情報を用いて状態変動正規分布モデルを更新してもよい。例えば、自動車等の移動体に本発明の位置算出装置を設定する場合において、車速検出システム（例えば車速パルス）や磁気センサー等の方位センサーから、移動体の速度や移動方向（つまり速度ベクトル）を参照情報として取得可能な場合は、当該参照情報を用いて状態変動正規分布モデルを更新するように構成することも可能である。
【０２２６】
具体的には、最新の移動体の位置ベクトル“Ｐ”に速度ベクトル“Ｖ”を加味して、現在の演算時刻における移動体の位置ベクトル“Ｐ”を算出する。このようにして位置ベクトル“Ｐ”が算出可能であるため、演算時刻毎の位置ベクトルの変化量“δＰ”を算出することができる。従って、速度ベクトル“Ｖ”を利用して位置変化量ベクトル“δＰ”を算出し、これを用いて状態変化量期待値ベクトル“μδX”を更新すればよい。
【０２２７】
また、状態変動正規分布モデルの標準偏差も同様に更新することが可能である。移動体の進行方向がわかれば、当該進行方向以外の方向に対する速度成分はほとんど無いため、当該方向に対する位置や速度の誤差はゼロに近くなると考えられる。そこで、状態変化量共分散行列“ＳδX”のうち、移動体の進行方向以外の方向に対応する共分散の成分を小さな値（例えばゼロ）とするように状態変化量共分散行列“ＳδX”を更新すると効果的である。
【０２２８】
５−４．非確率変数
上記の実施形態では、位置算出演算によって算出可能な算出パラメーターのうち、位置ベクトル“Ｐ”の各成分を確率変数とし、クロックバイアス“ｂ”を非確率変数とする場合を例に挙げて説明した。しかし、これはあくまでも一例であり、確率変数及び非確率変数とする算出パラメーターは適宜に選択可能である。
【０２２９】
例えば、３次元の位置成分のうちの一部の成分の変動を確率変数とするといったことも可能である。具体的には、位置変化量ベクトル“δＰ”のうちの高度方向の位置変化量“δｚ”を非確率変数とし、“δｚ”については受信信号を利用した所定の収束演算を行って近似解を求めることとしてもよい。また、高度方向の位置変化量“δｚ”及びクロックバイアス変化量“δｂ”の２種類の算出パラメーターを非確率変数とし、“δｚ”及び“δｂ”については収束演算を行って近似解を求めるといったことも当然に可能である。
【０２３０】
５−５．最尤推定アルゴリズム
上記の実施形態では、最尤推定のためのアルゴリズムとしてＥＭアルゴリズムを適用した場合を例に挙げて説明したが、他のアルゴリズムを適用してもよいことは勿論である。例えば、山登り法として知られる勾配法を適用してもよい。勾配法は、最適化問題に対するアルゴリズムの一種であり、所定の評価関数の勾配を用いてパラメーターの最適化を行う手法である。
【０２３１】
５−６．電子機器
上記の実施例では、ＧＰＳ受信機を具備する電子機器として携帯型電話機を例に挙げて説明したが、本発明を適用可能な電子機器はこれに限られるわけではない。例えば、カーナビゲーション装置や携帯型ナビゲーション装置、パソコン、ＰＤＡ（Personal Digital Assistant）、腕時計といった他の電子機器についても本発明を同様に適用可能である。
【０２３２】
５−７．衛星測位システム
また、上記の実施形態では、衛星測位システムとしてＧＰＳを例に挙げて説明したが、ＷＡＡＳ（Wide Area Augmentation System）、ＱＺＳＳ（Quasi Zenith Satellite System）、ＧＬＯＮＡＳＳ（GLObal NAvigation Satellite System）、ＧＡＬＩＬＥＯ等の他の衛星測位システムであってもよいことは勿論である。
【０２３３】
５−８．処理の主体
上記の実施例では、ＧＰＳ受信部の内部に設けられた処理部が位置算出処理を行うものとして説明した。つまり、上記の実施例では、ＧＰＳ受信部（ベースバンド処理回路部）が位置算出装置として機能するものとして説明した。しかし、位置算出の主体を、電子機器のプロセッサーであるホスト処理部としてもよい。
【０２３４】
この場合は、例えば、ベースバンド回路部の処理部がＧＰＳ衛星信号を捕捉してメジャメント情報を算出・取得し、ホスト処理部に出力する。そして、ホスト処理部は、処理部から入力したメジャメント情報を用いて、上記の原理に基づく位置算出演算を行う。この場合は、電子機器（ホスト処理部）が位置算出装置として機能することになる。
【符号の説明】
【０２３５】
１ 携帯型電話機、 ５ ＧＰＳアンテナ、 １０ ＧＰＳ受信部、 １１ ＲＦ受信回路部、 ２０ ベースバンド処理回路部、 ２１ 処理部、 ２３ 記憶部、 ３０ ホスト処理部、 ４０ 操作部、 ５０ 表示部、 ６０ 携帯電話用アンテナ、 ７０ 携帯電話用無線通信回路部、 ８０ 記憶部、 ９０ 時計部
【技術分野】
【０００１】
本発明は、位置算出方法及び位置算出装置に関する。
【背景技術】
【０００２】
測位用信号を利用した測位システムとしては、ＧＰＳ（Global Positioning System）が広く知られており、携帯型電話機やカーナビゲーション装置等に内蔵された位置算出装置に利用されている。ＧＰＳでは、複数のＧＰＳ衛星の位置や各ＧＰＳ衛星から位置算出装置までの擬似距離等の情報に基づいて、位置算出装置の位置座標及び時計誤差を求める位置算出演算を行う。
【０００３】
衛星測位システムを利用した位置算出演算では、複数の測位用衛星について受信信号を用いて測定した擬似距離に含まれ得る誤差（以下、「擬似距離誤差」と称す。）の二乗和を最小化させる、いわゆる最小二乗法を利用した位置算出演算が知られている（例えば、特許文献１。）。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００９−９７８９７号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
従来の位置算出演算では、位置算出演算を行う時刻（以下、「演算時刻」と称す。）毎に、最小二乗法を利用するなどして位置算出演算を１回１回個別に行うのが通常である。従って、例えば、測定した擬似距離に瞬間的に大きな誤差が混入した場合、そのときの位置算出演算の精度が急激に低下するという欠点があったが、以降の位置算出演算には影響しないという利点もあった。但し、瞬間的とは言え、真の位置から大きく離れた位置が演算結果として得られる、いわゆる位置飛びが発生することは、場合によって大きな問題となり得た。
【０００６】
本発明は上述した課題に鑑みて為されたものであり、その目的とするところは、位置算出の正確性を向上させるための新しい位置算出手法を提案することにある。
【課題を解決するための手段】
【０００７】
以上の課題を解決するための第１の形態は、測位用信号を受信することと、前記測位用信号の受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行うことと、を含む位置算出方法である。
【０００８】
また、他の形態として、測位用信号を受信する受信部と、前記受信部による受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う位置算出部と、を備えた位置算出装置を構成してもよい。
【０００９】
この第１の形態等によれば、測位用信号の受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う。これにより、各演算結果を独立とするのではなく、演算結果に時間的な相関を持たせて位置算出を行うことが可能となる。この場合、測定量に瞬間的に大きな誤差が混入したとしても、位置算出精度が低下することを防止できる。
【００１０】
また、第２の形態として、第１の形態の位置算出方法において、過去の前記位置算出演算の結果に基づいて前記変動の基準を更新することと、前記基準を用いて前記変動を算出することと、を更に含む、位置算出方法を構成することとしてもよい。
【００１１】
この第２の形態によれば、過去の位置算出演算の結果に基づいて演算結果の変動の基準を更新することで、過去の位置算出演算の結果を今回の位置算出演算の結果に反映させることが可能となる。更新した基準を用いて演算結果の変動を算出することで、位置算出の正確性が向上する。
【００１２】
また、第３の形態として、第１又は第２の形態の位置算出方法において、前記位置算出演算は、少なくとも前記演算結果の変動を確率変数とした正規分布モデルに基づいて定められた演算であり、過去の前記位置算出演算の結果に基づいて前記正規分布モデルを更新することを更に含む、位置算出方法を構成することとしてもよい。
【００１３】
この第３の形態によれば、少なくとも演算結果の変動を確率変数とした正規分布モデルに基づいて定められた位置算出演算を行う。つまり、演算結果の変動の分布が正規分布に従うと仮定する。その上で、正規分布モデルを固定的に設定するのではなく、過去の位置算出演算の結果に基づいて正規分布モデルを更新することで、時間変化に対応した適切な正規分布モデルに基づく位置算出が可能となる。
【００１４】
また、第４の形態として、第３の形態の位置算出方法において、前記受信信号に基づいて擬似距離を測定することを更に含み、前記位置算出演算は、前記演算結果の変動を確率変数とした正規分布モデルと、前記擬似距離に含まれる誤差を確率変数とした確率分布モデルとに基づいて定められた演算である、位置算出方法を構成することとしてもよい。
【００１５】
この第４の形態によれば、測位用信号の受信信号に基づいて擬似距離を測定する。そして、演算結果の変動を確率変数とした正規分布モデルと、擬似距離に含まれる誤差を確率変数とした確率分布モデルとに基づいて定められた位置算出演算を行う。演算結果の変動を確率変数とした正規分布モデルに加えて、擬似距離に含まれる誤差を確率変数とした確率分布モデルを用いることで、位置算出演算を適切に行うことができる。
【００１６】
より具体的には、第５の形態として、第４の形態の位置算出方法における前記擬似距離に含まれる誤差を確率変数とした確率分布モデルを、正規分布モデル又は正規混合分布モデルとする、位置算出方法を構成することができる。
【００１７】
通常の測位環境においては、擬似距離に含まれる誤差の分布は正規分布に従うと仮定することができる。しかし、例えばマルチパス環境のように、擬似距離に含まれる誤差の分布を正規分布と仮定することが必ずしも適切でない場合もある。そこで、第５の形態のように、擬似距離に含まれる誤差を確率変数とした確率分布モデルを、正規分布モデル又は正規混合分布モデルとして位置算出演算を行うことで、測位環境に見合った適切な位置算出を実現することが可能となる。
【００１８】
また、第６の形態として、第３〜第５の何れかの形態の位置算出方法において、前記位置算出演算は、位置情報の他に、クロックバイアス、速度情報及びクロックドリフトのうちの何れかを算出する演算であり、前記正規分布モデルは、当該演算結果の各要素の変動を確率変数としたモデルである、位置算出方法を構成することとしてもよい。
【００１９】
この第６の形態によれば、位置情報の他に、クロックバイアス、速度情報及びクロックドリフトのうちの何れかを算出する位置算出演算を行う。その際、演算結果の各要素の変動を確率変数とした正規分布モデルを用いることで、各要素に時間的な相関を持たせて位置算出を行うことが可能となる。
【００２０】
また、第７の形態として、第１〜第６の何れかの形態の位置算出方法において、前記位置算出演算によって算出可能な算出パラメーターのうちの所定のパラメーターの値を、前記受信信号を用いた所定の収束演算を行って算出することを更に含み、前記位置算出演算は、過去に前記収束演算によって算出された前記所定のパラメーターの値を用いる演算であることを含む、位置算出方法を構成することとしてもよい。
【００２１】
この第７の形態によれば、位置算出演算によって算出可能な算出パラメーターのうちの所定のパラメーターの値を、受信信号を用いた所定の収束演算を行って算出する。そして、過去に収束演算によって算出された所定のパラメーターの値を用いる位置算出演算を行う。これにより、特定の算出パラメーターについては、時間的な相関を持たせずに、所定の収束演算を行って値を算出することができる。そして、その算出結果を、所与の確率分布モデルに基づいて定められた位置算出演算に反映させることで、位置算出を効果的に行うことが可能となる。
【図面の簡単な説明】
【００２２】
【図１】状態変動正規分布モデルの説明図。
【図２】擬似距離誤差正規分布モデルの説明図。
【図３】マルチパスの説明図。
【図４】自己相関の一例を示す図。
【図５】マルチパス環境における自己相関の一例を示す図。
【図６】マルチパス環境における自己相関の一例を示す図。
【図７】第１の擬似距離誤差正規混合分布モデルの説明図。
【図８】第２の擬似距離誤差正規混合分布モデルの説明図。
【図９】位置算出結果の一例を示す図。
【図１０】位置算出結果の一例を示す図。
【図１１】位置算出結果の一例を示す図。
【図１２】携帯型電話機の機能構成の一例を示すブロック図。
【図１３】ベースバンド処理回路部の回路構成の一例を示す図。
【図１４】ベースバンド処理の流れを示すフローチャート。
【図１５】第１の位置算出処理の流れを示すフローチャート。
【図１６】第２の位置算出処理の流れを示すフローチャート。
【図１７】第３の位置算出処理の流れを示すフローチャート。
【図１８】第４の位置算出処理の流れを示すフローチャート。
【図１９】ドップラー測位処理の流れを示すフローチャート。
【発明を実施するための形態】
【００２３】
以下、図面を参照して、本発明の好適な実施形態の一例について説明する。本実施形態は、衛星測位システムの一種であるＧＰＳ（Global Positioning System）を適用した実施形態である。なお、本発明を適用可能な実施形態が以下説明する実施形態に限定されるわけでないことは勿論である。
【００２４】
１．原理
本実施形態における位置算出方法（位置算出の原理）について説明する。ＧＰＳを利用した衛星測位システムにおいて、測位用衛星の一種であるＧＰＳ衛星は、エフェメリスやアルマナックといった衛星軌道データを含む航法データを、測位用の衛星信号の一種であるＧＰＳ衛星信号に乗せて発信している。
【００２５】
ＧＰＳ衛星信号は、拡散符号の一種であるＣ／Ａ（Coarse and Acquisition）コードによって、スペクトラム拡散方式として知られるＣＤＭＡ（Code Division Multiple Access）方式によって変調された１．５７５４２［ＧＨｚ］の通信信号である。Ｃ／Ａコードは、コード長１０２３チップを１ＰＮフレームとする繰返し周期１ｍｓの擬似ランダム雑音符号であり、各ＧＰＳ衛星に固有のコードである。
【００２６】
ＧＰＳでは、複数のＧＰＳ衛星の位置や各ＧＰＳ衛星からＧＰＳ受信機（位置算出装置）までの擬似距離といった、ＧＰＳ衛星信号の受信信号に係る測定量を用いて、ＧＰＳ受信機の位置（位置座標）や時計誤差（クロックバイアス）を演算する位置算出演算を行う。擬似距離は、衛星サーチを行ってＧＰＳ衛星信号を捕捉し、その結果として得られるメジャメント情報を用いて測定される。
【００２７】
衛星サーチには、ＧＰＳ衛星信号を受信した信号に対する位相方向の相関演算（以下、「位相サーチ」と称す。）と、周波数方向の相関演算（以下、「周波数サーチ」と称す。）とが含まれる。位相サーチを行うことで、ＧＰＳ受信機がＧＰＳ衛星信号を受信した際のＣ／Ａコードの位相（以下、「コード位相」と称す。）をメジャメント情報として取得できる。また、周波数サーチを行うことで、ＧＰＳ受信機が受信したＧＰＳ衛星信号の受信周波数やドップラー周波数をメジャメント情報として取得できる。
【００２８】
観念的には、ＧＰＳ衛星とＧＰＳ受信機との間には複数のＣ／Ａコードが並んでいると考えることができる。ＧＰＳ衛星とＧＰＳ受信機との間の距離は、ちょうどＣ／Ａコードの整数倍の長さになるとは限らず、端数部分が生じ得る。この擬似距離の端数部分に相当するのがコード位相である。擬似距離の整数部分は、ＧＰＳ受信機及びＧＰＳ衛星の概略位置から求められる。
【００２９】
本実施形態では、ＧＰＳ受信機（位置算出装置）の位置ベクトル“Ｐ”及びクロックバイアス（時計誤差）“ｂ”を未知数とし、これらを成分とするＧＰＳ受信機の状態ベクトル“Ｘ”を次式（１）のように定義する。
【数１】
【００３０】
式（１）において、“（ｘ，ｙ，ｚ）”は、測位を行う座標系におけるＧＰＳ受信機（位置算出装置）の位置ベクトル“Ｐ”の各軸の位置成分である。後々の便宜のため、位置成分“（ｘ，ｙ，ｚ）”及びクロックバイアス“ｂ”を、状態“Ｘ”の成分という意味で“Ｘi＝（ｘ，ｙ，ｚ，ｂ）T＝（Ｘ1，Ｘ2，Ｘ3，Ｘ4）T”と表記する。但し、下付きの添え字の“ｉ＝１，２，３，４”は、状態“Ｘ”の成分の番号を示す。上付きの添え字の“Ｔ”は行列の転置を示す。また、本実施形態の数式では、ベクトルを意味する矢印の表記を省略して簡潔に記載する。
【００３１】
ＧＰＳ受信機の真の位置ベクトルを“Ｐ＝Ｐ0＋δＰ”、真のクロックバイアスを“ｂ＝ｂ0＋δｂ”とおく。但し、“δＰ”及び“δｂ”は、それぞれ位置ベクトル及びクロックバイアスの初期値に適用される未知の変化量（修正量）であり、それぞれ「位置変化量ベクトル」及び「クロックバイアス変化量」と呼称する。また、位置変化量ベクトル“δＰ”及びクロックバイアス変化量“δｂ”を纏めて、次式（２）のように状態変化量ベクトル“δＸ”を定義する。
【数２】
【００３２】
従来から用いられているＧＰＳの測位計算では、例えば次式（３）の観測方程式を利用して、状態変化量ベクトル“δＸ”を求める。
【数３】
【００３３】
本明細書では、ＧＰＳ衛星のうちＧＰＳ受信機が捕捉に成功した衛星（以下、「捕捉衛星」と称す。）の番号を“ｋ”と表記し、原則的に各諸量について下付きの添え字で表記する。また、本実施形態では、ＧＰＳ受信機が合計Ｋ個のＧＰＳ衛星の捕捉に成功したものとして説明する。つまり、「ｋ＝１，２，・・・，Ｋ」として説明する。
【００３４】
式（３）の左辺の“δρ”は擬似距離修正量ベクトルであり、次式（４）のように表される。
【数４】
つまり、擬似距離修正量ベクトル“δρ”は、各捕捉衛星の擬似距離修正量“δρk”を成分とするベクトルである。
【００３５】
また、式（３）の右辺の行列“Ｇ”は、捕捉衛星の衛星配置を決定付ける幾何行列であり、次式（５）ように表される。
【数５】
但し、“ｌk”は、ＧＰＳ受信機からｋ番目の捕捉衛星に向かう視線方向のベクトル（以下、「視線ベクトル」と称す。）である。
【００３６】
また、式（３）の右辺の“Ｅ”は擬似距離誤差ベクトルであり、次式（６）のように表される。
【数６】
但し、“εk”は、ｋ番目の捕捉衛星からＧＰＳ受信機までの擬似距離に含まれる誤差（以下、「擬似距離誤差」と称す。）である。擬似距離誤差ベクトル“Ｅ”は、各捕捉衛星の擬似距離誤差“εk”を成分とするベクトルである。
【００３７】
式（３）における未知数は、位置変化量ベクトル“δＰ＝（δｘ，δｙ，δｚ）”及びクロックバイアス変化量“δｂ”の合計４個である。従って、式（３）は、Ｋ≧４の場合に解くことができる。Ｋ＞４の場合は、式（３）は過決定方程式となる。この場合は、Ｋ個の捕捉衛星の擬似距離誤差“ε”の二乗和を最小とする解を求める最小二乗法が一般的に用いられる。具体的には、式（３）の最小二乗解は、例えば次式（７）のように求められる。
【数７】
【００３８】
式（３）の観測方程式を利用した位置算出演算では、演算時刻毎に、ＧＰＳ衛星信号を受信して擬似距離を測定し、式（７）の最小二乗解を解析的に求める。この場合、状態変化量ベクトル“δＸ”を、演算時刻毎に独立として計算している。
【００３９】
しかし、ＧＰＳ受信機（位置算出装置）は、通常は連続的に移動する。例えば、ＧＰＳ受信機を携行したユーザーが通常歩行や通常走行しているのであれば、演算時刻毎のＧＰＳ受信機の移動距離にそれほど大きな差は現れないはずである。つまり、ユーザーが突飛な動作をしない限り、ＧＰＳ受信機の位置変化量ベクトル“δＰ”の変動はほぼ一定となることが期待される。ＧＰＳ受信機を自動車等の移動体に設置した場合も同様である。このことから、本願発明者は、位置算出演算の演算結果の変動を確率変数とした所与の確率分布モデルを定義することを考えた。
【００４０】
例えば、ユーザーが通常歩行しており、「１秒間隔」で位置算出演算を行うとすると、１秒毎のＧＰＳ受信機の移動距離は「数メートル程度」となり、個人差を考慮しても、その誤差幅は「数十センチ程度」になると考えられる。従って、歩行を想定した場合、位置算出演算の演算結果のうち、位置変化量ベクトル“δＰ”については、期待値が「数メートル程度」、標準偏差が「数十センチ程度」と仮定できるであろう。また、ＧＰＳ受信機を自動車に設置する場合を考えると、位置変化量ベクトル“δＰ”は、期待値が「数十メートル程度」、標準偏差が「数メートル程度」と仮定できるであろう。
【００４１】
クロックバイアス“ｂ”についても同様に考えると、短い時間間隔で考えた場合、クロックバイアス“ｂ”がそれほど大きく変動するとは考えにくい。そのため、クロックバイアス変化量“δｂ”についても確率分布モデルを仮定できるであろう。
【００４２】
そこで、本願発明者は、位置ベクトル“Ｐ”及びクロックバイアス“ｂ”を成分とする状態ベクトル“Ｘ”について、その変動である状態変化量ベクトル“δＸ”を確率変数とする正規分布型の確率分布モデル（以下、「状態変動正規分布モデル」と称す。）を仮定した。
【００４３】
図１は、状態変動正規分布モデル（状態変動ＧＭ（Gaussian Model））の説明図である。図１には、状態変化量ベクトル“δＸ”の１つの成分“δＸi”に着目したモデルを図示している。横軸は状態変化量の成分“δＸi”であり、縦軸はその確率密度“ｐ（δＸi）”である。状態変化量“δＸi”の期待値“μ”を「状態変化量期待値」と呼称し、状態変化量“δＸi”の標準偏差“σ”を「状態変化量標準偏差」と呼称する。
【００４４】
状態変化量“δＸi”を確率変数とし、その分布が正規分布Ｎ（μ，σ2）に従うと仮定する。しかし、本実施形態では、状態変化量ベクトル“δＸ”には、位置変化量ベクトル“Ｐ”の各成分“（δｘ，δｙ，δｚ）”及びクロックバイアス変化量“δｂ”が含まれる。そこで、一次元の正規分布モデルを拡張して、次式（８）で与えられる多次元状態変動正規分布モデル“ｐδX(m)（δＸ）”を用いることとする。
【数８】
【００４５】
式（８）において、上付きの括弧書きの“ｍ”は確率変数として取り扱う算出パラメーターの数を示す。算出パラメーターとは、位置算出演算によって算出可能なパラメーターであり、本実施形態では“（δｘ，δｙ，δｚ，δｂ）”の４つのパラメーターがこれに相当する。４つの算出パラメーター全てを確率変数とするのであれば“ｍ＝４”となる。詳細は後述するが、４つの算出パラメーターのうちの一部の算出パラメーターを確率変数とすることも可能である。この場合は、確率変数とする算出パラメーターの数が“ｍ”となる。
【００４６】
式（８）の多次元状態変動正規分布モデル“ｐδX(m)（δＸ）”は、非独立のｍ個の確率変数の同時発生確率を与える確率密度関数である。“μδX”は、状態変化量ベクトル“δＸ”の各成分の期待値でなる状態変化量期待値ベクトルである。また、“ＳδX”は、状態変化量ベクトル“δＸ”の成分間の共分散でなる状態変化量共分散行列である。これらについては、数式を用いて詳細に後述する。
【００４７】
本実施形態では、上記の多次元状態変動正規分布モデル“ｐδX(m)（δＸ）”を適用して、最尤推定法を利用して状態変化量ベクトル“δＸ”の最尤推定解を求める。これは、測位用信号を受信し、測位用信号の受信信号を用いた演算であって、位置算出演算の演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行うことに相当する。状態変化量ベクトル“δＸ”の最尤推定解が求まれば、当該最尤推定解を状態ベクトルの初期値“Ｘ0”に加算することで、真の状態ベクトル“Ｘ”が推定できる。
【００４８】
さて、本実施形態の位置算出演算では、ＧＰＳ衛星信号の受信信号に基づいて擬似距離を測定する。そして、上記の状態変動正規分布モデルと、擬似距離に含まれる誤差を確率変数とした確率分布モデルとに基づいて定められた位置算出演算を行う。つまり、擬似距離誤差“ε”も確率変数として取り扱い、所与の確率分布モデルを仮定して演算を行う。本実施形態では、擬似距離誤差“ε”を確率変数とした確率分布モデルを、正規分布モデル又は正規混合分布モデルとする。
【００４９】
図２は、擬似距離誤差正規分布モデル（擬似距離誤差ＧＭ）の説明図である。図２において、横軸は擬似距離誤差“ε”であり、縦軸はその確率密度“ｐε（ε）”である。擬似距離誤差“ε”の期待値“με”を「擬似距離誤差期待値」と呼称し、擬似距離誤差“ε”の標準偏差“σε”を「擬似距離誤差標準偏差」と呼称する。
【００５０】
擬似距離誤差“ε”が、正規分布Ｎ（με，（σε）2）に従うと仮定して観測方程式の解を求める手法が最小二乗法である。広く用いられるのは「με＝０，σε＝１」とする擬似距離誤差正規分布モデルであり、この場合の最小二乗解が式（７）と一致する。
【００５１】
上記の擬似距離誤差正規分布モデル“ｐε（ε）”は、ＧＰＳ受信機の天空が開けた環境（例えばオープンスカイ環境）において有用である。オープンスカイ環境では、ＧＰＳ衛星から品質の良いＧＰＳ衛星信号を受信することができることから、擬似距離誤差“ε”はゼロ近傍の値となることが期待されるためである。しかし、ＧＰＳ受信機の天空が開けていない環境では、擬似距離誤差“ε”の挙動が変化する。その典型例はマルチパス環境である。
【００５２】
図３は、マルチパス環境の説明図である。直接波信号を点線で示し、間接波信号を一点鎖線で示す。間接波信号の存在により、ＧＰＳ受信機が測定する擬似距離には誤差が生ずる。擬似距離は、ＧＰＳ受信機が位相サーチを行って取得したコード位相を用いて算出されるが、マルチパス環境では、このコード位相に誤差が重畳する。コード位相の誤差は、真のコード位相に対して、正の誤差となる場合及び負の誤差となる場合がある。
【００５３】
図４は、コード位相検出の原理の説明図である。図４では、横軸をコード位相、縦軸を相関値として、Ｃ／Ａコードの自己相関値の概略例を示している。なお、以下の説明では、相関値というときは、相関値の大きさ（絶対値）を意味するものとする。
【００５４】
Ｃ／Ａコードの自己相関値は、理想的にはピーク値を頂点とする左右対称の略三角形の形状で表される。この場合に、相関値のピーク値（以下、「相関ピーク値」と称す。）に対応する位相が、受信したＧＰＳ衛星信号のＣ／Ａコードの位相である。相関ピーク値を検出するために、例えば、あるコード位相に対して、一定量だけ進んだコード位相（以下、「進み位相」と称す。）と、一定量だけ遅れたコード位相（以下、「遅れ位相」と称す。）とを考える。そして、進み位相の相関値と遅れ位相の相関値とが等しくなる相関値を相関ピーク値として検出する。
【００５５】
図４では、相関値は左右対称の略三角形の形状であるため、進み位相の相関値と遅れ位相の相関値との中心のコード位相の相関値が相関ピーク値となり、対応するコード位相が「ピーク位相」として検出される。この検出されたピーク位相のことを「検出ピーク位相」と称する。図４の相関値の形状は理想形状であるが、マルチパス環境では相関値の形状が変化する。
【００５６】
図５及び図６は、マルチパス環境における相関値の形状の一例を示す図である。図５は、間接波信号が直接波信号と同位相で到達した場合の相関値のグラフの一例であり、図６は、間接波信号が直接波信号と逆位相で到達した場合の相関値のグラフの一例である。これらの図には、直接波信号と、間接波信号と、この直接波信号と間接波信号とを合成したマルチパス信号とのそれぞれに対応する相関値のグラフを示している。横軸はコード位相、縦軸は相関値である。
【００５７】
間接波信号に対する相関値は、直接波信号に対する相関値と同様に略三角形の形状をなしているが、間接波信号の相関ピーク値の大きさは、直接波信号の相関ピーク値よりも小さい。これは、ＧＰＳ衛星から送出されたＧＰＳ衛星信号が、建物や地面に反射したり障害物を透過することによって、送出された時点における信号強度が、受信時には弱められていることによるものである。
【００５８】
また、間接波信号のピーク位相は、直接波信号のピーク位相よりも遅れている。これは、ＧＰＳ衛星から送出されたＧＰＳ衛星信号が、建物や地面に反射したり障害物を回折することによって、ＧＰＳ衛星からＧＰＳ受信機までの伝搬距離が長くなったことによるものである。マルチパス信号に対する相関値は、直接波信号の相関値と間接波信号の相関値との和となるため、三角形状が歪んでピーク値を中心とした左右対称とはならない。
【００５９】
間接波信号が直接波信号と同位相でＧＰＳ受信機に到達した場合は、直接波信号と間接波信号とは互いに強め合う。このため、合成波信号の相関値は、直接波信号に対する相関値の大きさと間接波信号に対する相関値の大きさとの合算値となる。この場合、相関値の形状は、例えば図５に示すような形状となる。図５において、検出ピーク位相は、真のピーク位相よりも遅れた位相となる。
【００６０】
それに対し、間接波信号が直接波信号よりも、例えば半周期以上１周期未満の範囲等の逆位相となった場合は、直接波信号と間接波信号とは互いに弱め合う。このため、合成波信号の相関値は、直接波信号に対する相関値の大きさから間接波信号に対する相関値の大きさを減じた値となる。この場合、相関値の形状は、例えば図６に示すような形状となる。図６において、検出ピーク位相は、真のピーク位相よりも進んだ位相となる。
【００６１】
ここで、検出ピーク位相と真のピーク位相との位相差を「コード位相誤差」と定義する。そして、便宜的に、検出ピーク位相が真のピーク位相よりも遅れている場合のコード位相誤差の符号を「正」、検出ピーク位相が真のピーク位相よりも進んでいる場合のコード位相誤差の符号を「負」と定義する。この場合、図５ではコード位相誤差は「正」となり、図６ではコード位相誤差は「負」となる。
【００６２】
前述したように、擬似距離はコード位相を用いて算出されるため、コード位相の誤差は擬似距離に誤差として重畳される。すなわち、図５のようにコード位相に正の誤差が含まれると、擬似距離誤差“ε”は正の誤差となる。一方、図６のようにコード位相に負の誤差が含まれると、擬似距離誤差“ε”は負の誤差となる。
【００６３】
このように、マルチパス環境等においては、コード位相誤差の正負の符号に応じて、擬似距離誤差“ε”は正の値にもなり得るし、負の値にもなり得る。勿論、擬似距離誤差“ε”にはコード位相誤差以外の誤差成分も含まれる。そのため、これらの誤差成分の大きさ如何によっても、擬似距離誤差“ε”は種々の値をとり得る。
【００６４】
このことから、本願発明者は、例えばマルチパス環境においては、擬似距離誤差“ε”の分布を正規分布と仮定するのではなく、複数の正規分布を混合した正規混合分布と仮定することが適切であると判断した。擬似距離誤差“ε”に正規混合分布を仮定した確率分布モデルとして、本実施形態では２種類のモデルを例示する。
【００６５】
図７は、擬似距離誤差ＧＭＭの一例である第１の擬似距離誤差正規混合分布モデル（第１の擬似距離誤差ＧＭＭ（Gaussian Mixture Model））の説明図である。第１の擬似距離誤差ＧＭＭ“ｆ（ε）”は、第１の正規分布モデル“ｐ1（ε）”と、第２の正規分布モデル“ｐ2（ε）”との２種類の正規分布モデルを合成した関数として定義される。図７において、横軸は擬似距離誤差“ε”（単位はメートル）であり、縦軸はその確率密度“ｆ（ε）”である。また、第１の正規分布モデル“ｐ1（ε）”を点線で示し、第２の正規分布モデル“ｐ2（ε）”を一点鎖線で示し、第１の擬似距離誤差ＧＭＭ“ｆ（ε）”を太実線でそれぞれ示す。
【００６６】
第１の正規分布モデル“ｐ1（ε）”は、ＧＰＳ受信機が直接波信号を受信することを想定した擬似距離誤差“ε”の確率分布モデルである。間接波信号の影響がなければ、擬似距離誤差“ε”は“０メートル”の近傍の値となることが期待される。そのため、擬似距離誤差期待値を“０メートル”（με＝０）とし、擬似距離誤差標準偏差を比較的小さな値である“１０メートル”（σε＝１０）とする正規分布関数“Ｎ（０，１０2）”によって、第１の正規分布モデル“ｐ1（ε）”をモデル化した。
【００６７】
第２の正規分布モデル“ｐ2（ε）”は、ＧＰＳ受信機が直接波信号及び間接波信号を受信することを想定した擬似距離誤差“ε”のモデル関数である。前述したように、コード位相誤差の正負の符号により、擬似距離誤差“ε”は正負の何れの値もとり得る。そのため、第２の正規分布モデル“ｐ2（ε）”の擬似距離誤差期待値を、正負の中心である“０メートル”（με＝０）とした。また、コード位相誤差の大きさ如何によっては、擬似距離誤差“ε”に非常に大きな誤差が含まれ得る。そのため、擬似距離誤差標準偏差を大きめに見積もって“１００メートル”（σε＝１００）とした。すなわち、第２の正規分布モデル“ｐ2（ε）”を正規分布関数“Ｎ（０，１００2）”によってモデル化した。
【００６８】
第１の正規分布モデル“ｐ1（ε）”と第２の正規分布モデル“ｐ2（ε）”とを混合する割合（合成比率）は、次のように決定した。ＧＰＳ受信機がマルチパスの影響を受ける割合は、一般的にはそれほど高くないと考えられる。そこで、ＧＰＳ受信機が１割の確率でマルチパスの影響を受けると仮定して、第１の正規分布モデル“ｐ1（ε）”の重みを“９／１０”（９割）とし、第２の正規分布モデル“ｐ1（ε）”の重みを“１／１０”（１割）とした。従って、第１の正規分布モデル“ｐ1（ε）”を９／１０倍した関数と、第２の正規分布モデル“ｐ2（ε）”を１／１０倍した関数とを混合（合成）した関数として、第１の擬似距離誤差ＧＭＭ“ｆ（ε）”をモデル化した。
【００６９】
図８は、擬似距離誤差ＧＭＭの別例である第２の擬似距離誤差正規混合分布モデル（第２の擬似距離誤差ＧＭＭ）の説明図である。第２の擬似距離誤差ＧＭＭ“ｆ（ε）”は、図７の第１の擬似距離誤差ＧＭＭと同様に、第１の正規分布モデル“ｐ1（ε）”と第２の正規分布モデル“ｐ2（ε）”との２種類の正規分布モデルを混合（合成）した確率密度関数である。
【００７０】
第１の正規分布モデル“ｐ1（ε）”は、ＧＰＳ受信機が直接波信号を受信することを想定した擬似距離誤差“ε”のモデル関数として、正規分布関数“Ｎ（０，１０2）”でモデル化した。また、第２の正規分布モデル“ｐ2（ε）”は、マルチパス環境においてＧＰＳ受信機が間接波信号を受信することを想定した擬似距離誤差“ε”のモデル関数として、正規分布関数“Ｎ（１００，２０2）”でモデル化した。
【００７１】
図７の第１の擬似距離誤差ＧＭＭ“ｆ（ε）”と異なるのは、第２の正規分布モデル“ｐ2（ε）”を、直接波信号及び間接波信号を受信することを想定したモデル関数として定義するのではなく、間接波信号単体を受信することを想定したモデル関数として定義した点にある。これは、例えば、ＧＰＳ受信機の近傍に高層ビル等の障害物が存在する環境において起こり得る。ＧＰＳ受信機の天空が完全に開けていないために、ＧＰＳ衛星から送出された信号は障害物に反射し、間接波信号としてＧＰＳ受信機に到達する。しかし、障害物に遮られて、直接波信号がＧＰＳ受信機に直接到達しない環境である。アーバンキャニオン環境が典型例である。
【００７２】
本願発明者が実験を行った結果、アーバンキャニオン環境では、擬似距離誤差“ε”として「＋１００〜＋１５０メートル」程度の正の誤差が含まれる傾向があることがわかった。そこで、第２の擬似距離誤差ＧＭＭ“ｆ（ε）”では、第２の正規分布モデル関数“ｐ2（ε）”の擬似距離誤差期待値を“１００メートル”（με＝１００）とし、擬似距離誤差標準偏差を“２０メートル”（σε＝２０）とした。すなわち、“ｐ2（ε）＝Ｎ（１００，２０2）”としてモデル化した。
【００７３】
第１及び第２の正規分布モデルの合成比率は、例えばＧＰＳ受信機が２割の確率でマルチパスの影響を受けると仮定して決定した。すなわち、第１の正規分布モデル“ｐ1（ε）”の重みを“８／１０”（８割）とし、第２の正規分布モデル“ｐ2（ε）”の重みを“２／１０”（２割）とした。従って、第２の正規分布モデル“ｐ1（ε）”を８／１０倍した関数と、第２の正規分布モデル“ｐ2（ε）”を２／１０倍した関数とを混合（合成）した関数として、第２の擬似距離誤差ＧＭＭ“ｆ（ε）”をモデル化した。
【００７４】
上記のように、擬似距離誤差“ε”を確率変数とする確率分布モデルとして、正規分布モデル及び正規混合分布モデルの２種類のモデルを定義した。位置算出を行う際には、測位環境に応じて、擬似距離誤差正規分布モデルと、擬似距離誤差正規混合分布モデルとを使い分けて位置算出演算を行う。以下、それぞれの確率分布モデルを適用する場合の具体的な位置算出方法について、数式を用いて詳細に説明する。
【００７５】
１−１．擬似距離誤差正規分布モデル（擬似距離誤差ＧＭ）を適用する場合
（１）第１の位置算出方法
第１の位置算出方法は、擬似距離誤差正規分布モデルを適用し、状態変化量ベクトル“δＸ”の全ての算出パラメーターを確率変数として位置算出する方法である。この方法では、状態変化量ベクトル“δＸ”の４つの算出パラメーター“（δｘ，δｙ，δｚ，δｂ）”全てを確率変数とする。この場合、式（８）の多次元状態変動正規分布モデル“ｐδX(m)（δＸ）”は、４次元状態変動正規分布モデル“ｐδX(4)（δＸ）”となる。
【００７６】
最尤推定法は、尤度と呼ばれる指標値に基づいて、所与の観測データを利用して推定対象とするパラメーターの最尤解を推定する統計学的手法である。尤度は、ある観測データについて、適用する確率モデルがどの程度当てはまっているかを確率で表す尺度である。
【００７７】
本実施形態において、尤度Ｌは、式（６）の擬似距離誤差ベクトル“Ｅ”のＫ衛星の擬似距離誤差“εk”（ｋ＝１，２，・・・，Ｋ）の同時発生確率と、状態変化量ベクトル“δＸ”の４つの算出パラメーター“δＸi＝（δｘ，δｙ，δｚ，δｂ）”（ｉ＝１，２，３，４）の同時発生確率との積で表される。つまり、尤度Ｌは、次式（９）のように表される。
【数９】
【００７８】
式（９）において、“ｐε,k（ε）”はｋ番目の捕捉衛星の擬似距離誤差正規分布モデルを示す。なお、図２で定義した擬似距離誤差正規分布モデル“ｐε（ε）”は、全ての捕捉衛星についてパラメーター値を共通化した共通のモデルとしてもよいし、捕捉衛星毎にパラメーター値を個別化した異なるモデルとしてもよい。
【００７９】
また、式（９）において、“ｐδX,i（δＸi）”はｉ番目の算出パラメーターの状態変動正規分布モデルを示す。１行目から２行目に変形する際に、４つの算出パラメーターの同時発生確率を与える４次元状態変動正規分布モデル“ｐδX(4)（δＸ）”に書き換えた。
【００８０】
式（９）の尤度Ｌを最大化する状態変化量ベクトル“δＸ”を最尤推定解として求める。尤度Ｌをそのまま用いて計算するのは困難であるため、本実施形態では、その対数をとった対数尤度ｌｏｇＬを用いて計算する。対数尤度ｌｏｇＬは、次式（１０）のように計算される。
【数１０】
【００８１】
式（１０）において、“ｃｏｎｓｔ”は定数項である。また、“μE”は、Ｋ衛星の擬似距離誤差“ε”の期待値“με,k＝（με,1，με,2，・・・，με,K）”を成分とする擬似距離誤差期待値ベクトルであり、次式（１１）で与えられる。
【数１１】
【００８２】
“ＷE”は、擬似距離誤差逆共分散行列であり、Ｋ衛星分の擬似距離誤差“εk”の共分散でなる擬似距離誤差共分散行列の逆行列として計算される。具体的には、次式（１２）で与えられる。
【数１２】
【００８３】
“μδX”は、式（８）で登場した状態変化量期待値ベクトルであり、次式（１３）で与えられる。
【数１３】
但し、“μδx，μδy，μδz，μδb”は、それぞれ“δｘ，δｙ，δｚ，δｂ”の期待値である。
【００８４】
また、“ＳδX”は、式（８）で登場した状態変化量共分散行列であり、次式（１４）で与えられる。
【数１４】
但し、４×４の各成分は、“δｘ，δｙ，δｚ，δｂ”の４成分間の共分散である。
【００８５】
また、“ＷδX”は、状態変化量逆共分散行列であり、状態変化量共分散行列“ＳδX”の逆行列として計算される。具体的には、次式（１５）で与えられる。
【数１５】
【００８６】
式（１０）を最大化する“δＸ”を求めるためには、“ｄ（ｌｏｇＬ）／ｄ（δＸ）＝０”を解けばよい。途中の計算は省略するが、最終的に、最尤推定解“δＸ”は次式（１６）のように求められる。
【数１６】
但し、“δρ”は式（４）の擬似距離変化量ベクトルであり、“Ｇ”は式（５）の幾何行列である。
【００８７】
各演算時刻それぞれについて、式（１６）に従って状態変化量ベクトル“δＸ”の最尤推定解を求める。そして、求めた状態変化量ベクトル“δＸ”の最尤推定解を用いて、状態ベクトル“Ｘ”を算出する。
【００８８】
また、本実施形態の位置算出演算では、上記の最尤推定解の算出と併せて、過去の位置算出演算の結果に基づいて状態変動正規分布モデルを更新する処理を行う。ここで、演算時刻“ｔ”における状態変化量ベクトル“δＸ”を“δＸ（ｔ）”と表記する。また、演算時刻“ｔ”における状態変化量期待値ベクトル“μδX”及び状態変化量共分散行列“ＳδX”を、それぞれ“μδX（ｔ）”及び“ＳδX（ｔ）”と表記する。
【００８９】
この場合、状態変化量期待値ベクトル“μδX（ｔ）”及び状態変化量共分散行列“ＳδX（ｔ）”を、それぞれ次式（１７）及び（１８）に従って更新する。
【数１７】
【数１８】
【００９０】
式（１７）は、１つ前の演算時刻に計算された状態変化量ベクトル“δＸ（ｔ）”で、今回の演算時刻における状態変化量期待値ベクトル“μδX（ｔ＋１）”を更新することを意味する。これは、過去の位置算出演算の結果に基づいて、演算結果（状態変化量ベクトル“δＸ”）の変動の基準（状態変化量期待値ベクトル“μδX（ｔ））を更新することに相当する。また、式（１８）は、１つ前の演算時刻に計算された状態変化量共分散行列“ＳδX（ｔ）”等の諸量を用いて、今回の演算時刻における状態変化量共分散行列“ＳδX（ｔ＋１）”を更新することを意味する。
【００９１】
状態変化量期待値ベクトル“μδX”は、状態変動正規分布モデルの期待値に係るベクトルである。また、状態変化量共分散行列“ＳδX”は、状態変動正規分布モデルの標準偏差（共分散）に係る行列である。よって、式（１７）及び（１８）は、過去の位置算出演算の結果に基づいて状態変動正規分布モデルを更新する式であると言える。
【００９２】
（２）第２の位置算出方法
第２の位置算出方法は、擬似距離誤差“ε”に正規分布モデルを適用し、状態変化量ベクトル“δＸ”の一部の算出パラメーターを確率変数として位置算出する方法である。この方法では、状態変化量ベクトル“δＸ”を構成する４つの算出パラメーター“δｘ，δｙ，δｚ，δｂ”のうちの一部の算出パラメーターを確率変数とし、残余の算出パラメーターを非確率変数として取り扱う。
【００９３】
確率変数及び非確率変数とする算出パラメーターは、適用するシステムに応じて適宜選択可能である。ここでは、ユーザーの位置変化量ベクトル“δＰ”を構成する“（δｘ，δｙ，δｚ）”の３つの算出パラメーターを確率変数とし、クロックバイアス変化量“δｂ”を非確率変数とする場合を例に挙げて説明する。この場合、式（８）の多次元状態変動正規分布モデル“ｐδX(m)（δＸ）”は、３次元状態変動正規分布モデル“ｐδX(3)（δＸ）”となる。
【００９４】
この場合、クロックバイアス変化量“δｂ”については正規分布モデルを仮定しないため、式（１０）に相当する対数尤度“ｌｏｇＬ”は、次式（１９）のように変更される。
【数１９】
【００９５】
また、式（１３）で定義した状態変化量期待値ベクトル“μδX”は、次式（２０）のように変更される。
【数２０】
但し、クロックバイアス変化量“δｂ”を非確率変数とするため、それに対応する４番目の成分をゼロとしている。
【００９６】
また、式（１４）で定義した状態変化量共分散行列“ＳδX”は、次式（２１）のように変更される。
【数２１】
但し、クロックバイアス変化量“δｂ”を非確率変数とするため、４×４の行列の４行目及び４列目の成分を全てゼロとしている。
【００９７】
この場合も、最尤推定解“δＸ”を式（１６）のように求めることができる。また、式（２０）で定義した状態変化量期待値ベクトル“μδX”と、式（２１）で定義した状態変化量共分散行列“ＳδX”とを、それぞれ式（１７）及び式（１８）に従って更新する。この場合、クロックバイアス変化量“δｂ”に対応する成分はゼロとして更新しない。
【００９８】
また、非確率変数とするクロックバイアス変化量“δｂ”については、受信信号を用いた所定の収束演算を行って解を求める。具体的には、式（３）の観測方程式に対して擬似距離誤差正規分布モデルを適用して、例えばニュートン法を利用した収束演算を行って、クロックバイアス変化量“δｂ”の近似解を求めることになる。収束演算によって算出されたクロックバイアス変化量“δｂ”は、次回の演算時刻での位置算出演算において擬似距離を近似的に求めるために使用する。
【００９９】
１−２．擬似距離誤差正規混合分布モデル（擬似距離誤差ＧＭＭ）を適用する場合
（３）第３の位置算出方法
第３の位置算出方法は、擬似距離誤差“ε”に正規混合分布モデルを適用し、状態変化量ベクトル“δＸ”の全ての算出パラメーターを確率変数として位置算出する方法である。考え方は第１の位置算出方法と同じであるが、擬似距離誤差の確率分布モデルとして、正規分布モデルではなく、正規混合分布モデルを適用する点が異なる。
【０１００】
この場合、擬似距離誤差ＧＭＭ“ｆ（ε）”と、４次元状態変動正規分布モデル“ｐδX(4)（δＸ）”とを用いて、対数尤度“ｌｏｇＬ”を次式（２２）のように定式化する。
【数２２】
【０１０１】
式（２２）において、“ｆk（εk）”は、ｋ番目の捕捉衛星の擬似距離誤差ＧＭＭを示し、次式（２３）のように一般化される。
【数２３】
【０１０２】
式（２３）において、“ｌ”は、混合（合成）する正規分布モデルの番号であり、Ｌ個（ｌ＝１，２，・・・，Ｌ）の正規分布モデルを混合する場合の式として記載している。“αl”は、ｌ番目の正規分布モデル“Ｎl”の重み（合成割合）である。“μk(l)”は、ｋ番目の捕捉衛星に係るｌ番目の正規分布モデルの期待値であり、“σk(l)”は、ｋ番目の捕捉衛星に係るｌ番目の正規分布モデルの標準偏差である。
【０１０３】
なお、上記の擬似距離誤差ＧＭＭは“ｆ（ε）”は、全ての捕捉衛星についてパラメーター値を共通化した共通のモデルとしてもよいし、捕捉衛星毎にパラメーター値を個別化した異なるモデルとしてもよい。
【０１０４】
式（２２）の対数尤度“ｌｏｇＬ”を最大化する“δＸ”を求めるためには、“ｄ（ｌｏｇＬ）／ｄ（δＸ）＝０”を解けばよい。しかし、正規混合分布（ＧＭＭ）では、対数の加算など複雑な項が存在するため、解析的に解を求めることが困難である。そこで、本実施形態では、最適化問題に対するアルゴリズムの一種であるＥＭ（Expectation Maximum）アルゴリズムを用いて最尤推定解を求める。
【０１０５】
ＥＭアルゴリズムは、確率分布モデルのパラメーターを最尤推定法に基づいて推定する手法の一種である。観測可能なデータ（サンプル）の他に、観測不可能なデータ（隠しパラメーター）が存在する場合であって、当該隠しパラメーターが確率分布モデルに依存する場合に用いられる手法である。ＥＭアルゴリズムは反復法の一種であり、尤度関数の条件付き確率に関する期待値を計算するＥ（Expectation）ステップと、期待値を最大化する解を求めるＭ（Maximum）ステップとの２つのステップを反復して最尤推定解を求める。
【０１０６】
途中の計算は省略するが、ＥＭアルゴリズムを適用すると、状態変化量ベクトル“δＸ”の最尤推定解は、次式（２４）のように求まる。
【数２４】
但し、“δＸ0”は各演算時刻における状態変化量ベクトル“δＸ”の初期値である。
【０１０７】
式（２４）における“Ｍl（δＸ0）”は、次式（２５）で与えられる。
【数２５】
【０１０８】
また、式（２５）における“ｈk（ｌ）（δＸ0）”は、次式（２６）で与えられる。
【数２６】
【０１０９】
また、式（２４）における“μl”は、次式（２７）で与えられる。
【数２７】
【０１１０】
各演算時刻において、式（２４）〜式（２７）に従って最尤推定解“δＸ”を求める。また、擬似距離誤差正規分布モデルを適用した場合と同様に、式（１３）で定義した状態変化量期待値ベクトル“μδX”と、式（１４）で定義した状態変化量共分散行列“ＳδX”とを、演算時刻毎に更新する。
【０１１１】
具体的には、状態変化量期待値ベクトル“μδX（ｔ）”及び状態変化量共分散行列“ＳδX（ｔ）”を、それぞれ次式（２８）及び（２９）に従って更新する。
【数２８】
【数２９】
【０１１２】
式（２９）における“Ｈ（δＸ）”及び“Ｆ（δＸ）”は、それぞれ次式（３０）及び（３１）で与えられる。
【数３０】
【数３１】
【０１１３】
また、式（３１）における“Ｒl”は、次式（３２）で与えられる。
【数３２】
【０１１４】
（４）第４の位置算出方法
第４の位置算出方法は、擬似距離誤差“ε”に正規混合分布モデルを適用し、状態変化量ベクトル“δＸ”の一部の算出パラメーターを確率変数として位置算出する方法である。この方法では、第２の位置算出方法と同様に、状態変化量ベクトル“δＸ”を構成する４つの算出パラメーター“δｘ，δｙ，δｚ，δｂ”のうちの一部の算出パラメーターを確率変数とし、残余の算出パラメーターを非確率変数として取り扱う。擬似距離誤差の確率分布モデルとして、正規分布モデルではなく、正規混合分布モデルを適用する点が第２の位置算出方法と異なる。
【０１１５】
この場合も、第３の位置算出方法と同様の数式を用いて、状態変化量ベクトル“δＸ”の最尤推定解を求めることができる。第３の位置算出方法と異なるのは、確率変数として取り扱う算出パラメーターについてのみ正規分布モデルを適用して計算する点である。つまり、第２の位置算出方法と同様に、状態変化量ベクトル“δＸ”の成分のうち、非確率変数として取り扱う算出パラメーターに対応する成分をゼロとして計算する。
【０１１６】
例えば、クロックバイアス変化量“δｂ”を非確率変数とする場合は、クロックバイアス変化量“δｂ”に対応する成分をゼロとして計算する。そして、クロックバイアス変化量“δｂ”については、受信信号を用いた所定の収束演算を行って解を算出する。具体的には、式（３）の観測方程式に対して擬似距離誤差正規混合分布モデル（擬似距離誤差ＧＭＭ）を適用して、例えばＥＭアルゴリズムを利用した収束演算を行って、クロックバイアス変化量“δｂ”の近似解を求めることになる。この場合、収束演算によって算出されたクロックバイアス変化量“δｂ”は、次回の演算時刻での位置算出演算において擬似距離を近似的に求めるために使用する。
【０１１７】
２．実験結果
本実施形態の位置算出方法で位置算出演算を行った実験結果について説明する。従来の手法と本実施形態の手法とのそれぞれを用いて、位置算出演算を行って実際に位置を算出する実験を行った。ユーザーに位置算出装置を携行させ、予め定められた周回経路に沿って歩行させながら位置を算出した。
【０１１８】
図９は、位置算出結果の一例を示す実験結果である。グラフの横方向は東西方向を示し、縦方向は南北方向を示す。単位はそれぞれメートル［ｍ］である。ここで示す結果は、状態変動の確率分布モデルとして正規分布を適用し、擬似距離誤差の確率分布モデルとして正規混合分布（ＧＭＭ）を適用した結果である。従来の位置算出方法で得られた軌跡を白丸のプロットで示し、本実施形態の位置算出方法で得られた軌跡を黒丸のプロットで示す。
【０１１９】
この結果を見ると、従来の位置算出方法を適用した場合は、所々に位置飛びが発生しており、それほど正確な位置は得られていないことがわかる。しかし、本実施形態の位置算出方法を適用した場合は、周回経路に沿った滑らかな軌跡が得られていることがわかる。位置算出の演算結果を独立とするのではなく、演算結果に時間的な相関を持たせて位置算出を行ったことで、算出位置の連続性が飛躍的に向上している。
【０１２０】
図１０は、上記の歩行試験において位置演算精度を計測したグラフである。グラフの横方向は東西方向の位置誤差を示し、縦方向は南北方向の位置誤差を示す。単位はそれぞれメートル［ｍ］である。グラフの中心が位置誤差“０［メートル］”に相当し、プロットがグラフの中心に近いほど位置算出精度が高いことを意味する。
【０１２１】
このグラフを見ると、従来の位置算出方法を適用した場合は、プロットが全体的に広くばらついており、位置算出精度が低くなっていることがわかる。それに対し、本実施形態の位置算出方法を適用した場合は、プロットが全体的に中心部に集中しており、従来の手法と比べて位置算出精度が格段に向上していることがわかる。
【０１２２】
図１１は、位置算出精度を示す他の実験結果の一例である。ここで示す結果は、状態変動の確率分布モデル及び擬似距離誤差の確率分布モデルとして、共に正規分布モデルを適用した結果である。図の見方は図１０と同じである。この結果を見ると、擬似距離誤差の確率分布モデルとして正規分布モデルを適用した場合も、位置誤差を示すプロットがグラフの中心部に集中しており、位置算出精度が格段に向上していることがわかる。
【０１２３】
上記の実験結果により、本実施形態の位置算出方法を用いることで、位置算出の正確性が向上することが実証された。
【０１２４】
３．実施例
次に、上記の位置算出方法を用いて位置算出を行う位置算出装置の実施例について説明する。ここでは、位置算出装置を具備した電子機器の一例である携帯型電話機の実施例として説明する。なお、本発明を適用可能な実施例が以下説明する実施例に限定されるわけでないことは勿論である。
【０１２５】
３−１．携帯型電話機の構成
図１２は、携帯型電話機１の機能構成の一例を示すブロック図である。携帯型電話機１は、ＧＰＳアンテナ５と、ＧＰＳ受信部１０と、ホスト処理部３０と、操作部４０と、表示部５０と、携帯電話用アンテナ６０と、携帯電話用無線通信回路部７０と、記憶部８０と、時計部９０とを備えて構成される。
【０１２６】
ＧＰＳアンテナ５は、ＧＰＳ衛星から発信されているＧＰＳ衛星信号を含むＲＦ（Radio Frequency）信号を受信するアンテナであり、受信信号をＧＰＳ受信部１０に出力する。
【０１２７】
ＧＰＳ受信部１０は、ＧＰＳアンテナ５から出力された信号に基づいて携帯型電話機１の位置を算出する回路或いは装置であり、ＧＰＳ受信機に相当する機能ブロックである。本実施例では、このＧＰＳ受信部１０が位置算出装置に相当する。
【０１２８】
ＧＰＳ受信部１０は、ＲＦ受信回路部１１と、ベースバンド処理回路部２０とを備えて構成される。ＧＰＳ受信部１０は、ＧＰＳ衛星信号を受信する受信部に相当する。なお、ＲＦ受信回路部１１と、ベースバンド処理回路部２０とは、それぞれ別のＬＳＩ（Large Scale Integration）として製造することも、１チップとして製造することも可能である。
【０１２９】
ＲＦ受信回路部１１は、ＲＦ信号の受信回路である。回路構成としては、例えば、ＧＰＳアンテナ５から出力されたＲＦ信号をＡ／Ｄ変換器でデジタル信号に変換し、デジタル信号を処理する受信回路を構成してもよい。また、ＧＰＳアンテナ５から出力されたＲＦ信号をアナログ信号のまま信号処理し、最終的にＡ／Ｄ変換することでデジタル信号をベースバンド処理回路部２０に出力する構成としてもよい。
【０１３０】
後者の場合には、例えば、次のようにＲＦ受信回路部１１を構成することができる。すなわち、所定の発振信号を分周或いは逓倍することで、ＲＦ信号乗算用の発振信号を生成する。そして、生成した発振信号を、ＧＰＳアンテナ５から出力されたＲＦ信号に乗算することで、ＲＦ信号を中間周波数の信号（以下、「ＩＦ（Intermediate Frequency）信号」と称す。）にダウンコンバートし、ＩＦ信号を増幅等した後、Ａ／Ｄ変換器でデジタル信号に変換して、ベースバンド処理回路部２０に出力する。
【０１３１】
ベースバンド処理回路部２０は、ＲＦ受信回路部１１から出力された受信信号に対して、キャリア除去や相関演算等を行ってＧＰＳ衛星信号を捕捉する。そして、捕捉したＧＰＳ衛星信号から抽出した時刻情報や衛星軌道情報等を利用して、携帯型電話機１の位置及び時計誤差を算出する。
【０１３２】
ホスト処理部３０は、記憶部８０に記憶されているシステムプログラム等の各種プログラムに従って携帯型電話機１の各部を統括的に制御するプロセッサーであり、ＣＰＵ（Central Processing Unit）等のプロセッサーを有して構成される。ホスト処理部３０は、ベースバンド処理回路部２０から取得した位置座標をもとに、表示部５０に現在位置を指し示した地図を表示させたり、その位置座標を各種のアプリケーション処理に利用する。
【０１３３】
操作部４０は、例えばタッチパネルやボタンスイッチ等により構成される入力装置であり、押下されたキーやボタンの信号をホスト処理部３０に出力する。この操作部４０の操作により、通話要求やメール送受信要求、位置算出要求等の各種指示入力がなされる。
【０１３４】
表示部５０は、ＬＣＤ（Liquid Crystal Display）等により構成され、ホスト処理部３０から入力される表示信号に基づいた各種表示を行う表示装置である。表示部５０には、位置表示画面や時刻情報等が表示される。
【０１３５】
携帯電話用アンテナ６０は、携帯型電話機１の通信サービス事業者が設置した無線基地局との間で携帯電話用無線信号の送受信を行うアンテナである。
【０１３６】
携帯電話用無線通信回路部７０は、ＲＦ変換回路、ベースバンド処理回路等によって構成される携帯電話の通信回路部であり、携帯電話用無線信号の変調・復調等を行うことで、通話やメールの送受信等を実現する。
【０１３７】
記憶部８０は、ＲＯＭ（Read Only Memory）やフラッシュＲＯＭ、ＲＡＭ（Random Access Memory）等の記憶装置を有して構成され、ホスト処理部３０が携帯型電話機１を制御するためのシステムプログラムや、各種アプリケーション処理を実行するための各種プログラムやデータ等を記憶する。
【０１３８】
時計部９０は、携帯型電話機１の内部時計であり、水晶振動子及び発振回路でなる水晶発振器等を有して構成される。時計部９０の計時時刻は、ベースバンド処理回路部２０及びホスト処理部３０に随時出力される。時計部９０は、ベースバンド処理回路部２０により算出された時計誤差を用いて較正される。
【０１３９】
３−２．ベースバンド処理回路部の構成
図１３は、ベースバンド処理回路部２０の回路構成及び記憶部２３のデータ構成の説明図である。ベースバンド処理回路部２０は、主要な機能構成として、処理部２１と、記憶部２３とを備える。
【０１４０】
処理部２１は、ベースバンド処理回路部２０の各機能部を統括的に制御する制御装置及び演算装置であり、ＣＰＵやＤＳＰ（Digital Signal Processor）等のプロセッサーを有して構成される。処理部２１は、衛星捕捉部２１１と、位置算出部２１３とを機能部として有する。
【０１４１】
衛星捕捉部２１１は、ＧＰＳ衛星の捕捉を行う機能部である。具体的には、ＲＦ受信回路部１１から出力されるデジタル化された受信信号に対して、キャリア除去及び相関演算のデジタル信号処理を実行して、ＧＰＳ衛星を捕捉する。そして、相関演算結果に対するピーク判定を行い、受信キャリア信号のドップラー周波数や受信周波数、受信Ｃ／Ａコードのコード位相といった情報をメジャメント情報として取得する。
【０１４２】
また、衛星捕捉部２１１は、相関演算結果に基づいて航法データを復号する。受信キャリア信号の位相（キャリア位相）と、受信Ｃ／Ａコードの位相（コード位相）とが検出され、相関がとれた状態になると、相関値の時間変化をもとに、航法データを構成する各ビット値をデコードすることができる。この位相同期は、例えば位相ロックループとして知られるＰＬＬ（Phase Locked Loop）処理により実現される。
【０１４３】
位置算出部２１３は、衛星捕捉部２１１により捕捉された各ＧＰＳ衛星に係るメジャメント情報や航法データ、時刻情報、衛星情報といった諸量を用いて、上記の原理に従って位置算出演算を行って、携帯型電話機１の位置（位置座標）及び時計誤差（クロックバイアス）を算出する。そして、算出した位置をホスト処理部３０に出力するとともに、算出した時計誤差で時計部９０を較正する。
【０１４４】
位置算出部２１３は、ＧＰＳ受信部１０による受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う位置算出部に相当する。
【０１４５】
記憶部２３は、ベースバンド処理回路部２０のシステムプログラムや、衛星捕捉機能、位置算出機能といった各種機能を実現するための各種プログラム、データ等を記憶する。また、各種処理の処理中データ、処理結果などを一時的に記憶するワークエリアを有する。
【０１４６】
記憶部２３には、プログラムとして、処理部２１により読み出され、ベースバンド処理（図１４参照）として実行される第１のベースバンド処理プログラム２３１が記憶されている。また、ベースバンド処理プログラム２３１は、各種の位置算出処理（図１５〜図１８参照）として実行される位置算出プログラム２３１１をサブルーチンとして含む。原理で説明した第１〜第４の位置算出方法の何れかに対応するプログラムが位置算出プログラム２３１１として記憶されている。
【０１４７】
また、記憶部２３には、主要なデータとして、衛星軌道データ２３２と、捕捉した各ＧＰＳ衛星に係る衛星別データ２３３と、状態変動正規分布モデル２３４と、演算結果データ２３５とが記憶される。
【０１４８】
衛星軌道データ２３２は、全てのＧＰＳ衛星の概略の衛星軌道情報を記憶したアルマナックや、各ＧＰＳ衛星それぞれについて詳細な衛星軌道情報を記憶したエフェメリス等のデータである。この衛星軌道データ２３２は、ＧＰＳ衛星から受信したＧＰＳ衛星信号をデコードすることで取得する他、例えば携帯型電話機１の基地局やアシストサーバーからアシストデータとして取得する。
【０１４９】
衛星別データ２３３は、メジャメント情報２３３Ａと、衛星情報２３３Ｂと、測位計算用諸量２３３Ｃとを含む。メジャメント情報２３３Ａは、当該捕捉衛星について測定したＧＰＳ衛星信号のコード位相や受信周波数、ドップラー周波数といった情報である。衛星情報２３３Ｂは、当該捕捉衛星の位置や速度、移動方向といった情報である。測位計算用諸量２３３Ｃは、上記の原理に従って測位計算を行う際に用いる各種の諸量である。
【０１５０】
状態変動正規分布モデル２３４は、図１で図示・説明した状態変動正規分布モデルであり、モデルパラメーター値として、例えば、状態変化量期待値２３４Ａと、状態変化量標準偏差２３４Ｂとを含む。状態変動正規分布モデルは、演算時刻毎に更新される。
【０１５１】
演算結果データ２３５は、位置算出演算の結果として得られたデータであり、携帯型電話機１の位置ベクトル“Ｐ”やクロックバイアス“ｂ”がこれに含まれる。
【０１５２】
３−３．処理の流れ
図１４は、記憶部２３に記憶されているベースバンド処理プログラム２３１が処理部２１により読み出されることで、ベースバンド処理回路部２０において実行されるベースバンド処理の流れを示すフローチャートである。
【０１５３】
最初に、処理部２１は、記憶部２３にサブルーチンとして記憶されている位置算出プログラム２３１１に従って位置算出処理を開始する（ステップＡ１）。位置算出処理は、図１５〜図１８の第１〜第４の位置算出処理の何れかの処理である。
【０１５４】
次いで、処理部２１は、演算結果の出力タイミングであるか否かを判定する（ステップＡ３）。演算結果の出力タイミングは、任意のタイミングを設定可能である。例えば、位置算出演算の時間間隔と同じ時間間隔のタイミングとしてもよいし、位置算出演算の時間間隔よりも長い時間間隔のタイミングとしてもよい。また、ユーザーから演算結果の出力が指示されたタイミングを出力タイミングとしてもよい。
【０１５５】
出力タイミングであると判定した場合は（ステップＡ３；Ｙｅｓ）、処理部２１は、最新の演算結果をホスト処理部３０に出力する（ステップＡ５）。ステップＡ３において出力タイミングではないと判定した場合（ステップＡ３；Ｎｏ）、又は、ステップＡ５の後、処理部２１は、ベースバンド処理を終了するか否かを判定する（ステップＡ７）。
【０１５６】
ベースバンド処理を終了しないと判定した場合は（ステップＡ７；Ｎｏ）、処理部２１は、ステップＡ３に戻る。また、終了すると判定した場合は（ステップＡ７；Ｙｅｓ）、処理部２１は、ベースバンド処理を終了する。
【０１５７】
図１５は、第１の位置算出方法に対応する第１の位置算出処理の流れを示すフローチャートである。最初に、処理部２１は、状態ベクトル“Ｘ”の初期値を設定する（ステップＢ１）。具体的には、サーバーアシストによって携帯型電話機１の基地局から取得した位置（基地局位置）や、予め定められた位置（例えば固定値）を、初期位置として設定する。また、クロックバイアス“ｂ”にも所定の初期値（例えば固定値）を設定する。
【０１５８】
次いで、処理部２１は、捕捉対象衛星選定処理を行う（ステップＢ３）。具体的には、時計部９０で計時されている現在時刻において、初期位置の天空に位置するＧＰＳ衛星を、記憶部２３に記憶されたアルマナックやエフェメリス等の衛星軌道データ２３２を用いて判定して、捕捉対象衛星に選定する。
【０１５９】
次いで、処理部２１は、状態変化量推定処理を行う（ステップＢ５〜Ｂ１７）。具体的には、処理部２１は、各捕捉対象衛星についてループＡを実行する（ステップＢ５〜Ｂ１３）。ループＡでは、処理部２１は、当該捕捉対象衛星からＧＰＳ衛星信号を受信して擬似距離を測定し、その値を擬似距離測定値“ρc”とする（ステップＢ７）。擬似距離測定値“ρc”は、位相サーチの結果として得られるコード位相を用いて算出することができる。
【０１６０】
また、処理部２１は、擬似距離近似値“ρa”を算出する（ステップＢ９）。具体的には、衛星軌道データ２３２から求まる当該捕捉対象衛星の位置ベクトル“Ｐk”と、携帯型電話機１の最新の位置ベクトル“Ｐ”及び最新のクロックバイアス“ｂ”とを用いて、擬似距離近似値を“ρa＝‖Ｐk−Ｐ‖＋ｂ”として算出する。但し、この式では、クロックバイアス“ｂ”の単位を距離に換算している。
【０１６１】
その後、処理部２１は、擬似距離変化量“δρ”を算出する（ステップＢ１１）。具体的には、ステップＢ７で測定した擬似距離測定値“ρc”と、ステップＢ９で算出した擬似距離近似値“ρa”との差を算出して、擬似距離変化量“δρ”とする。そして、処理部２１は、次の捕捉対象衛星へと処理を移行する。全ての捕捉対象衛星についてステップＢ７〜Ｂ１１を行った後、処理部２１は、ループＡを終了する（ステップＢ１３）。
【０１６２】
次いで、処理部２１は、幾何行列Ｇを算出する（ステップＢ１５）。具体的には、衛星軌道データ２３２から求まる各捕捉対象衛星の衛星位置ベクトルと、携帯型電話機１の最新の位置ベクトルとを用いて、携帯型電話機１から各捕捉対象衛星に向かう視線ベクトル“ｌ”を計算する。そして、式（５）の幾何行列Ｇを求める。
【０１６３】
その後、処理部２１は、擬似距離誤差正規分布モデル及び状態変動正規分布モデルを適用して、状態変化量“δＸ”の最尤推定解を求める（ステップＢ１７）。具体的には、式（１６）に従って状態変化量“δＸ”を計算する。擬似距離誤差“ε”に正規分布を仮定するため、ここで得られる最尤推定解は最小二乗解に相当する。そして、処理部２１は、状態変化量推定処理を終了する。
【０１６４】
次いで、処理部２１は、ステップＢ１７で最尤推定解として推定した状態変化量ベクトル“δＸ”を用いて状態ベクトル“Ｘ”を修正する（ステップＢ１９）。具体的には、最尤推定解として求めた状態変化量ベクトル“δＸ”を最新の状態ベクトル“Ｘ”に加算することで、今回の演算時刻における状態ベクトル“Ｘ”を新たに算出する。
【０１６５】
その後、処理部２１は、位置算出を終了するか否かを判定し（ステップＢ２３）、まだ終了しないと判定した場合は（ステップＢ２３；Ｎｏ）、状態変動正規分布モデル更新処理を行う（ステップＢ２５）。具体的には、状態変化量期待値ベクトル“μδX”及び状態変化量共分散行列“ＳδX”を、それぞれ式（１７）及び式（１８）に従って更新する。
【０１６６】
状態変動正規分布モデル更新処理を行った後、処理部２１は、ステップＢ３に戻り、次の演算時刻での位置算出演算へと移行する。また、ステップＢ２３において位置算出を終了すると判定した場合は（ステップＢ２３；Ｙｅｓ）、処理部２１は、第１の位置算出処理を終了する。
【０１６７】
ステップＢ３〜Ｂ２５までの一連の処理が、一の演算時刻での位置算出演算に相当する。第１の位置算出処理において特徴的であるのは、各演算時刻での位置算出演算において収束演算が省略されている点である。従来の位置算出方法では、最小二乗法を利用して解を求める際に、例えばニュートン法を利用した収束演算を行って近似解を求めていた。しかし、本実施形態では、演算結果の変動を確率変数とし、その上、演算時刻毎に確率分布モデルを更新することにしているため、収束演算を不要とする。
【０１６８】
図１６は、第２の位置算出方法に対応する第２の位置算出処理の流れを示すフローチャートである。状態変化量ベクトル“δＸ”のうち、位置変化量ベクトル“δＰ”の成分“δｘ，δｙ，δｚ”を確率変数とし、クロックバイアス変化量“δｂ”を非確率変数として位置算出を行うフローチャートである。なお、図１５の第１の位置算出処理と同一のステップについては同一の符号を付して説明を省略する。
【０１６９】
処理部２１は、ステップＢ３において捕捉対象衛星選定処理を行った後、状態変化量推定処理を行う（ステップＣ５）。具体的には、処理部２１は、第１の位置算出処理の状態変化量推定処理（ステップＢ５〜Ｂ１７）と同様の処理を行う。但し、ここではクロックバイアス“δｂ”を非確率変数とするため、状態変化量ベクトル“δＸ”のクロックバイアス変化量“δｂ”に対応する成分をゼロとして計算する。つまり、状態変化量推定処理は、位置変化量ベクトル“δＰ”を推定する位置変化量推定処理となる。
【０１７０】
また、状態変化量推定処理では、１回前の位置算出演算で収束演算によって算出したクロックバイアス“ｂ”の値を用いて、擬似距離近似値“ρa”を算出する（ステップＢ９）。具体的には、衛星軌道データ２３２から求まる当該捕捉対象衛星の位置ベクトル“Ｐk”と、前回の位置算出演算で算出した位置ベクトル“Ｐ”と、前回の位置算出演算で収束演算によって算出したクロックバイアス“ｂ”とを用いて、擬似距離近似値“ρa＝‖Ｐk−Ｐ‖＋ｂ”を算出する。これは、過去に収束演算によって算出された所定のパラメーターの値を用いて位置算出演算を行うことに相当する。
【０１７１】
次いで、処理部２１は、クロックバイアス変化量推定処理を行う（ステップＣ７〜Ｃ１９）。クロックバイアス変化量推定処理は、ＧＰＳ衛星信号の受信信号を用いた所定の収束演算を行って、非確率変数とするクロックバイアス変化量“δｂ”の近似解を求める処理である。
【０１７２】
詳細には、処理部２１は、各捕捉対象衛星についてループＢを行う（ステップＣ７〜Ｃ１３）。ループＢでは、処理部２１は、擬似距離近似値“ρa”を算出する（ステップＣ９）。そして、処理部２１は、擬似距離変化量“δρ”を算出する（ステップＣ１１）。具体的には、状態変化量推定処理のステップＢ７で測定した擬似距離測定値“ρc”と、ステップＣ９で算出した擬似距離近似値“ρa”との差を算出して、擬似距離変化量“δρ”とする。
【０１７３】
その後、処理部２１は、次の捕捉対象衛星へと処理を移行する。全ての捕捉対象衛星についてステップＣ９及びＣ１１を行った後、処理部２１は、ループＢを終了する（ステップＣ１３）。
【０１７４】
次いで、処理部２１は、擬似距離誤差正規分布モデルを適用して、クロックバイアス変化量“δｂ”の最小二乗解を求める（ステップＣ１７）。そして、処理部２１は、求めたクロックバイアス変化量“δｂ”について、所定の収束条件が成立したか否かを判定する（ステップＣ１８）。具体的には、例えば、今回の反復ステップで求めたクロックバイアス変化量“δｂ”と、１回前の反復ステップで求めたクロックバイアス変化量“δｂ”との差の絶対値が所定の閾値よりも小さくなった場合に、収束条件が成立したと判定する。
【０１７５】
収束条件が成立していないと判定した場合は（ステップＣ１８；Ｎｏ）、処理部２１は、ステップＣ１７で求めたクロックバイアス変化量“δｂ”を用いて、クロックバイアス“ｂ”を修正して更新する（ステップＣ１９）。そして、処理部２１は、ステップＣ７に戻る。また、収束条件が成立したと判定した場合は（ステップＣ１８；Ｙｅｓ）、処理部２１は、クロックバイアス変化量推定処理を終了する。
【０１７６】
その後、処理部２１は、状態変化量推定処理において推定した位置変化量ベクトル“δＰ”と、クロックバイアス変化量推定処理において推定したクロックバイアス変化量“δｂ”とを状態変化量ベクトル“δＸ”とし、これを用いて状態ベクトル“Ｘ”を修正する（ステップＣ２１）。つまり、最新の状態ベクトル“Ｘ”に状態変化量ベクトル“δＸ”を加算して、今回の演算時刻における状態ベクトル“Ｘ”を新たに算出する。
【０１７７】
その後、処理部２１は、位置算出を終了するか否かを判定し（ステップＢ２３）、まだ終了しないと判定した場合は（ステップＢ２３；Ｎｏ）、状態変動正規分布モデル更新処理を行う（ステップＢ２５）。その後、処理部２１は、ステップＢ３に戻り、次の演算時刻での位置算出演算へと移行する。また、ステップＢ２３において位置算出を終了すると判定した場合は（ステップＢ２３；Ｙｅｓ）、処理部２１は、第２の位置算出処理を終了する。
【０１７８】
この場合も、ステップＢ３〜Ｂ２５までの一連の処理が、一の演算時刻での位置算出演算に相当する。また、ステップＣ７〜Ｃ１９の繰り返し処理が、受信信号を用いた所定の収束演算に相当する。第２の位置算出処理において特徴的であるのは、確率変数として取り扱う位置変化量ベクトル“δＰ”については収束演算を省略しているが、非確率変数として取り扱うクロックバイアス変化量“δｂ”については収束演算を省略していない点である。クロックバイアス変化量“δｂ”は確率変数とせず、全くの未知数として取り扱うため、収束演算を行わなければ最適解が求まらないためである。
【０１７９】
図１７は、第３の位置算出方法に対応する第３の位置算出処理の流れを示すフローチャートである。第３の位置算出処理の基本的な流れは、図１５の第１の位置算出処理と同じである。しかし、第１の位置算出処理のステップＢ１７がステップＤ１７に置き換わる点が大きく異なる。
【０１８０】
ステップＤ１７では、処理部２１は、擬似距離誤差正規混合分布モデル及び状態変動正規分布モデルを適用して、状態変化量ベクトル“δＸ”の最尤推定解を求める。具体的には、処理部２１は、式（２４）〜式（２７）に従って、状態変化量ベクトル“δＸ”の最尤推定解を計算する。このステップは、ＥＭアルゴリズムを適用した計算ステップである。
【０１８１】
また、第３の位置算出処理では、ステップＢ２５の状態変動正規分布モデル更新処理において、処理部２１は、状態変化量期待値ベクトル“μδX”を式（２８）に従って更新する。また、処理部２１は、状態変化量共分散行列“ＳδX”を式（２９）〜式（３２）に従って更新する。
【０１８２】
図１８は、第４の位置算出方法に対応する第４の位置算出処理の流れを示すフローチャートである。第４の位置算出処理の基本的な流れは、図１６の第２の位置算出処理と同じである。しかし、第２の位置算出処理のステップＣ５がステップＥ５に置き換わり、ステップＣ１７がステップＥ１７に置き換わる点が大きく異なる。
【０１８３】
第４の位置算出処理では、ステップＥ５の状態変化量推定処理において、クロックバイアス“δｂ”を非確率変数として取り扱い、状態変化量ベクトル“δＸ”のクロックバイアス変化量“δｂ”に対応する成分をゼロとして計算を行う。つまり、この場合の状態変化量推定処理は、位置変化量ベクトル“δＰ”を推定する位置変化量推定処理となる。
【０１８４】
また、第２の位置算出処理と同様に、状態変化量推定処理では、１回前の位置算出演算で収束演算によって算出されたクロックバイアス“ｂ”の値を用いて、擬似距離近似値“ρa”を算出する（ステップＢ９）。つまり、過去に収束演算によって算出された所定のパラメーターの値を用いて位置算出演算を行う。
【０１８５】
また、第４の位置算出処理では、ステップＣ７〜ステップＣ１９のクロックバイアス変化量推定処理において、擬似距離誤差正規混合分布モデルを適用して、例えばＥＭアルゴリズムを利用した収束演算を行って（ステップＣ１８；Ｎｏ→ステップＣ１９→ステップＣ７）、クロックバイアス変化量“δｂ”の最尤推定解を求める（ステップＥ１７）。
【０１８６】
そして、処理部２１は、クロックバイアス変化量“δｂ”について所定の収束条件が成立したか否かを判定し（ステップＣ１８）、成立していないと判定した場合は（ステップＣ１８；Ｎｏ）、ステップＥ１７で求めたクロックバイアス変化量“δｂ”を用いてクロックバイアス“ｂ”を修正・更新した後（ステップＣ１９）、ステップＣ７に戻る。また、収束条件が成立したと判定した場合は（ステップＣ１８；Ｙｅｓ）、処理部２１は、ステップＣ２１へと移行する。
【０１８７】
４．作用効果
本実施形態によれば、測位用衛星信号の一種であるＧＰＳ衛星信号を受信する。そして、ＧＰＳ衛星信号の受信信号を用いた演算であって、ＧＰＳ受信機の位置及びクロックバイアスの変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う。
【０１８８】
具体的には、ＧＰＳ受信機（位置算出装置）の位置及びクロックバイアスの変動を確率変数とする状態変動正規分布モデルを設定する。また、擬似距離に含まれる誤差の変動を確率変数とする確率分布モデルとして、擬似距離誤差正規分布モデル又は擬似距離誤差正規混合分布モデルを設定する。そして、これらの確率モデルを用いて、最尤推定法に基づく演算を行って、位置及びクロックバイアスの変化量の最尤推定解を求める。これにより、演算時刻毎の演算結果を独立とするのではなく、演算結果に時間的な相関を持たせて位置算出を行うことが可能となる。この場合、測定した擬似距離に瞬間的に大きな誤差が混入したとしても、位置算出精度が低下することを防止できる。
【０１８９】
また、本実施形態では、過去の位置算出演算の結果に基づいて演算結果の変動の基準を更新し、更新した基準を用いて演算結果の変動を算出する。詳細には、１つ前の演算時刻に算出された位置及びクロックバイアスの変化量を用いて、今回の演算時刻での位置及びクロックバイアスの変化量の基準を更新する。これは、状態変動正規分布モデルの期待値を更新することに相当する。
【０１９０】
また、本実施形態では、１つ前の演算時刻での位置及びクロックバイアスの変化量の共分散を用いて、今回の演算時刻での位置及びクロックバイアスの変化量の共分散を更新する。これは、状態変動正規分布モデルの標準偏差を更新することに相当する。つまり、演算結果の変動に係る正規分布モデルを固定的に設定するのではなく、１つ前の演算時刻の演算結果に基づいて当該正規分布モデルを更新する。これにより、位置及びクロックバイアスの時間変化に俊敏に対応することができるようになり、移動体の位置をより正確に求めることが可能となる。
【０１９１】
５．変形例
５−１．ドップラー測位
上記の実施形態では、擬似距離を用いた位置算出方法を例に挙げて説明したが、ドップラーを利用した位置算出方法に対しても、本発明を実質的に同一に適用可能である。ドップラーを利用した位置算出方法は「ドップラー測位」と呼ばれる。
【０１９２】
ドップラーに関して、次式（３３）のモデル式が成立する。
【数３３】
但し、単位を全て速度の次元［ｍ／ｓ］として式を記載している。
【０１９３】
式（３３）において、“ｈk”は、ｋ番目の衛星について測定した周波数のずれ（以下、「周波数ずれ」と称す。）である。つまり、ＧＰＳ衛星信号の受信信号の周波数（受信周波数）と規定された搬送波周波数（１．５７５４２［ＧＨｚ］）との差として計算される。受信周波数は、周波数方向の相関演算を行うことで取得することができる。
【０１９４】
“Ｖk”は、ｋ番目の捕捉衛星の速度ベクトルであり、“Ｖ”は、ＧＰＳ受信機の速度ベクトルである。従って、“ＶR＝Ｖk−Ｖ”は、ＧＰＳ受信機とｋ番目の捕捉衛星との相対速度ベクトルである。また、“ｌk”は、ＧＰＳ受信機からｋ番目の捕捉衛星に向かう視線ベクトルである。相対速度ベクトル“ＶR＝Ｖk−Ｖ”を視線ベクトル“ｌk”の方向に投影することで、視線方向の相対速度ベクトル（ドップラーシフト）が求まる。
【０１９５】
また、式（３３）において、“ｄ”は、ＧＰＳ受信機の動作クロックのドリフト（クロックドリフト）である。また、“εkh”は、ｋ番目の衛星の周波数ずれに内在する観測誤差である。
【０１９６】
式（３３）に対してニュートン法を適用することを考える。ＧＰＳ受信機の位置ベクトル“Ｐ”、クロックバイアス“ｂ”、速度ベクトル“Ｖ”及びクロックドリフト“ｄ”を成分とする状態ベクトル“Ｘ＝（Ｐ，ｂ，Ｖ，ｄ）T”を定義する。つまり、位置情報の他に、クロックバイアス、速度情報及びクロックドリフトを状態ベクトル“Ｘ”の成分とする。
【０１９７】
ここで、式（３３）の（右辺）−（左辺）で表される次式（３４）の関数“ｇ（Ｐ，Ｖ）”を定義する。
【数３４】
【０１９８】
関数“ｇ”は、速度ベクトル“Ｖ”及び視線ベクトル“ｌ”を含む。視線ベクトル“ｌ”は、位置ベクトル“Ｐ”を変数とする。そのため、関数“ｇ”は、位置ベクトル“Ｐ”及び速度ベクトル“Ｖ”を変数とする関数である。
【０１９９】
関数“ｇ（Ｐ，Ｖ）”の位置ベクトル“Ｐ”、速度ベクトル“Ｖ”、クロックバイアス“ｂ”及びクロックドリフト“ｄ”での偏微分を求めると、それぞれ次式（３５）〜（３８）のようになる。
【数３５】
【数３６】
【数３７】
【数３８】
【０２００】
ここで、捕捉したＫ個の衛星について関数“ｇ（Ｐ，Ｖ）”を纏めた関数“Ｇ（Ｐ，Ｖ）”を便宜的に観測関数と呼び、次式（３９）のように定義する。
【数３９】
【０２０１】
この場合、次式（４０）の観測方程式が成立する。
【数４０】
【０２０２】
式（４０）の左辺は、位置ベクトルの初期値“Ｐ0”と、クロックバイアスの初期値“ｂ0”と、速度ベクトルの初期値“Ｖ0”と、クロックドリフトの初期値“ｄ0”とにより定まる観測関数の値“Ｇ（Ｐ0，ｂ0，Ｖ0，ｄ0）”である。
【０２０３】
また、式（４０）の右辺において、ベクトル“（δＰ，δｂ，δＶ，δｄ）T”を状態ベクトル“Ｘ”の変動を示す状態変化量ベクトル“δＸ”と定義する。また、便宜的に、観測関数“Ｇ”の各変数での偏微分で定義される行列を偏微分行列“Ｙ”と定義する。状態変化量ベクトル“δＸ”には、各変数の初期値により定まる偏微分行列“Ｙ”の値が乗算される。
【０２０４】
このように、式（３３）から、状態変化量ベクトル“δＸ＝（δＰ，δｂ，δＶ，δｄ）T”を未知数とする式（４０）の観測方程式が導出される。この観測方程式を利用して、上記の実施形態と同様に、状態変化量ベクトル“δＸ”の最尤推定解を求める。上記の実施形態と異なるのは、位置算出演算の基礎となる観測方程式が、式（３）から式（４０）に置き換わる点である。
【０２０５】
この場合、状態変化量ベクトル“δＸ＝（δＰ，δｂ，δＶ，δｄ）T”に含まれる算出パラメーター“（δｘ，δｙ，δｚ，δｂ，δｖx，δｖy，δｖz，δｄ）”の全部又は一部を確率変数とする確率分布モデルを定義する。位置情報の他に、クロックバイアス、速度情報及びクロックドリフトの変動を確率変数として取り扱う。
【０２０６】
８つの算出パラメーターの全てを確率変数とする場合は、式（８）の多次元状態変動正規分布モデルを、８次元状態変動正規分布モデル“ｐδX(8)（δＸ）”として定義する。また、８つの算出パラメーターのうち、例えばクロックバイアス“ｂ”及びクロックドリフト“ｄ”を除いた６つの算出パラメーターを確率変数とするのであれば、式（８）の多次元状態変動正規分布モデルを、６次元状態変動正規分布モデル“ｐδX(6)（δＸ）”として定義する。
【０２０７】
なお、式（４０）では未知数が８個であるため、解を求めるためには、少なくとも８個の衛星を捕捉する必要がある。そこで、未知数を減らす工夫をすると実用的である。例えば、速度ベクトル“Ｖ”を固定値（例えばゼロ）とし、位置ベクトル“Ｐ”、クロックバイアス“ｂ”及びクロックドリフト“ｄ”を未知数とする。この場合、式（４０）の観測方程式は、次式（４１）のように書き換えることができる。
【数４１】
【０２０８】
式（４１）では、状態ベクトル“Ｘ”の変動を示す状態変化量ベクトルを“δＸ＝（δＰ，δｂ，δｄ）T”と定義している。速度ベクトル“Ｖ＝（ｖx，ｖy，ｖz）”を固定値とすることで、未知数を８個から５個まで減らすことができる。これにより、少なくとも５個の衛星を捕捉することができれば、式（４１）の解を求めることができる。
【０２０９】
図１９は、上記の実施例における処理部２１が実行するドップラー測位処理の流れを示すフローチャートである。このドップラー測位処理は、図１４のベースバンド処理のステップＡ１において位置算出処理に置き換えて適用可能な処理である。また、このドップラー測位処理は、式（４０）の観測方程式に基づく処理であり、状態変化量ベクトル“δＸ”の全ての算出パラメーターを確率変数とする図１５の第１の位置算出処理に対応する処理である。
【０２１０】
最初に、処理部２１は、状態ベクトル“Ｘ”の初期値を設定する（ステップＦ１）。つまり、位置ベクトル“Ｐ”、クロックバイアス“ｂ”、速度ベクトル“Ｖ”及びクロックドリフト“ｄ”に所定の初期値を設定する。
【０２１１】
次いで、処理部２１は、捕捉対象衛星選定処理を行う（ステップＦ３）。そして、処理部２１は、状態変化量推定処理を行う（ステップＦ５〜Ｆ１９）。具体的には、処理部２１は、各捕捉対象衛星についてループＣを実行する（ステップＦ５〜Ｆ１３）。
【０２１２】
ループＣでは、処理部２１は、当該捕捉対象衛星からＧＰＳ衛星信号を受信して周波数ずれ“ｈ”を算出する（ステップＦ７）。周波数ずれ“ｈ”は、周波数サーチを行うことで取得した受信信号の周波数（受信周波数）と、規定された搬送波の周波数（規定搬送波周波数）との差として算出される。この周波数ずれ“ｈ”は、ドップラー測位における観測量となる。
【０２１３】
次いで、処理部２１は、位置算出装置と当該捕捉対象衛星との相対速度ベクトル“ＶR”を算出する（ステップＦ９）。相対速度ベクトル“ＶR”は、位置算出装置の最新の速度ベクトル“Ｖ”と、当該捕捉対象衛星の最新の速度ベクトル“Ｖk”とを用いて算出される。初回測位の場合は、速度ベクトルの初期値“Ｖ0”を最新の速度ベクトルとし、２回目以降の測位の場合は、前回の測位で算出した速度ベクトルを最新の速度ベクトルとして計算する。
【０２１４】
その後、処理部２１は、位置算出装置から当該捕捉対象衛星に向かう視線ベクトル“ｌ”を算出する（ステップＦ１１）。視線ベクトル“ｌ”は、位置算出装置の最新の位置ベクトル“Ｐ”と、当該捕捉対象衛星の最新の位置ベクトル“Ｐk”とを用いて算出される。初回測位の場合は、位置ベクトルの初期値“Ｐ0”を最新の位置ベクトルとし、２回目以降の測位の場合は、前回の測位で算出した位置ベクトルを最新の位置ベクトルとして計算する。
【０２１５】
その後、処理部２１は、次の捕捉対象衛星へと処理を移行する。全ての捕捉対象衛星についてステップＦ７〜Ｆ１１を行った後、処理部２１は、ループＣを終了する（ステップＦ１３）。
【０２１６】
ループＣの後、処理部２１は、最新の位置ベクトル“Ｐ”、クロックバイアス“ｂ”、速度ベクトル“Ｖ”及びクロックドリフト“ｄ”を用いて、式（４０）の左辺の観測関数“Ｇ”の値を算出する（ステップＦ１５）。また、処理部２１は、最新の位置ベクトル“Ｐ”、クロックバイアス“ｂ”、速度ベクトル“Ｖ”及びクロックドリフト“ｄ”を用いて、式（４０）の右辺の偏微分行列“Ｙ”の値を算出する（ステップＦ１７）。
【０２１７】
次いで、処理部２１は、式（４０）の観測方程式に従って、観測誤差正規分布モデル及び状態変動正規分布モデルを適用して、状態変化量“δＸ”の最尤推定解を求める（ステップＦ１９）。観測誤差正規分布モデルは、観測量である周波数ずれ“ｈ”に含まれる観測誤差“εh”を確率変数とし、その分布として正規分布を仮定したモデルである。
【０２１８】
ステップＦ１９を行った後、処理部２１は、状態変化量推定処理を終了する。以降のステップは、第１の位置算出処理と同様である。
【０２１９】
なお、ここではフローチャートの図示及び説明を省略するが、図１６〜図１８の第２〜第４の位置算出処理に対応するドップラー測位処理についても、それぞれ同様に構成することが可能である。
【０２２０】
５−２．確率変数
上記の実施形態では、位置算出演算の演算結果の変動を確率変数とする確率分布モデルを定義したが、位置算出演算の演算結果そのものを確率変数とする確率分布モデルを定義して位置算出演算を行うことも可能である。
【０２２１】
位置ベクトル“Ｐ”に着目して考えてみる。位置ベクトル“Ｐ”は、位置ベクトルの初期値“Ｐ0”に位置変化量ベクトル“δＰ”を加算することで求められる。位置ベクトルの初期値“Ｐ0”は固定値であるため、位置変化量ベクトル“δＰ＝Ｐ−Ｐ0”を確率変数として取り扱うことと、位置ベクトル“Ｐ”そのものを確率変数として取り扱うこととは同値である。従って、位置変化量ベクトル“δＰ”の代わりに位置ベクトル“Ｐ”そのものを確率変数として位置算出演算を行うことも可能である。クロックバイアス“ｂ”、速度ベクトル“Ｖ”及びクロックドリフト“ｄ”についても同様である。
【０２２２】
この場合は、例えばＧＰＳ受信機の状態ベクトル“Ｘ”そのものを確率変数とする状態正規分布モデルを定義する。そして、当該状態正規分布モデルを用いて、上記の実施形態と同様に位置算出演算を行う。状態正規分布モデルの更新では、過去に算出された位置ベクトル“Ｐ”やクロックバイアス“ｂ”を用いて、状態正規分布モデルの期待値や標準偏差を更新する。
【０２２３】
５−３．状態変動正規分布モデルの更新
上記の実施形態では、１時刻前の位置算出演算の結果に基づいて状態変動正規分布モデルを更新したが、あくまでもこれは一例である。例えば、２時刻前や３時刻前の位置算出演算の結果に基づいて状態変動正規分布モデルを更新してもよい。
【０２２４】
また、過去所定期間分（例えば５時刻分や１０時刻分）の位置算出演算の結果を平均処理し、その結果に基づいて状態変動正規分布モデルを更新してもよい。平均処理としては、単純な算術平均や幾何平均を適用してもよいし、加重平均を適用してもよい。加重平均を適用する場合は、現在の演算時刻により近い時刻での位置算出演算の結果ほど、重みをより高く設定して加重平均すると効果的である。
【０２２５】
また、外部から取得可能な参照情報を用いて状態変動正規分布モデルを更新してもよい。例えば、自動車等の移動体に本発明の位置算出装置を設定する場合において、車速検出システム（例えば車速パルス）や磁気センサー等の方位センサーから、移動体の速度や移動方向（つまり速度ベクトル）を参照情報として取得可能な場合は、当該参照情報を用いて状態変動正規分布モデルを更新するように構成することも可能である。
【０２２６】
具体的には、最新の移動体の位置ベクトル“Ｐ”に速度ベクトル“Ｖ”を加味して、現在の演算時刻における移動体の位置ベクトル“Ｐ”を算出する。このようにして位置ベクトル“Ｐ”が算出可能であるため、演算時刻毎の位置ベクトルの変化量“δＰ”を算出することができる。従って、速度ベクトル“Ｖ”を利用して位置変化量ベクトル“δＰ”を算出し、これを用いて状態変化量期待値ベクトル“μδX”を更新すればよい。
【０２２７】
また、状態変動正規分布モデルの標準偏差も同様に更新することが可能である。移動体の進行方向がわかれば、当該進行方向以外の方向に対する速度成分はほとんど無いため、当該方向に対する位置や速度の誤差はゼロに近くなると考えられる。そこで、状態変化量共分散行列“ＳδX”のうち、移動体の進行方向以外の方向に対応する共分散の成分を小さな値（例えばゼロ）とするように状態変化量共分散行列“ＳδX”を更新すると効果的である。
【０２２８】
５−４．非確率変数
上記の実施形態では、位置算出演算によって算出可能な算出パラメーターのうち、位置ベクトル“Ｐ”の各成分を確率変数とし、クロックバイアス“ｂ”を非確率変数とする場合を例に挙げて説明した。しかし、これはあくまでも一例であり、確率変数及び非確率変数とする算出パラメーターは適宜に選択可能である。
【０２２９】
例えば、３次元の位置成分のうちの一部の成分の変動を確率変数とするといったことも可能である。具体的には、位置変化量ベクトル“δＰ”のうちの高度方向の位置変化量“δｚ”を非確率変数とし、“δｚ”については受信信号を利用した所定の収束演算を行って近似解を求めることとしてもよい。また、高度方向の位置変化量“δｚ”及びクロックバイアス変化量“δｂ”の２種類の算出パラメーターを非確率変数とし、“δｚ”及び“δｂ”については収束演算を行って近似解を求めるといったことも当然に可能である。
【０２３０】
５−５．最尤推定アルゴリズム
上記の実施形態では、最尤推定のためのアルゴリズムとしてＥＭアルゴリズムを適用した場合を例に挙げて説明したが、他のアルゴリズムを適用してもよいことは勿論である。例えば、山登り法として知られる勾配法を適用してもよい。勾配法は、最適化問題に対するアルゴリズムの一種であり、所定の評価関数の勾配を用いてパラメーターの最適化を行う手法である。
【０２３１】
５−６．電子機器
上記の実施例では、ＧＰＳ受信機を具備する電子機器として携帯型電話機を例に挙げて説明したが、本発明を適用可能な電子機器はこれに限られるわけではない。例えば、カーナビゲーション装置や携帯型ナビゲーション装置、パソコン、ＰＤＡ（Personal Digital Assistant）、腕時計といった他の電子機器についても本発明を同様に適用可能である。
【０２３２】
５−７．衛星測位システム
また、上記の実施形態では、衛星測位システムとしてＧＰＳを例に挙げて説明したが、ＷＡＡＳ（Wide Area Augmentation System）、ＱＺＳＳ（Quasi Zenith Satellite System）、ＧＬＯＮＡＳＳ（GLObal NAvigation Satellite System）、ＧＡＬＩＬＥＯ等の他の衛星測位システムであってもよいことは勿論である。
【０２３３】
５−８．処理の主体
上記の実施例では、ＧＰＳ受信部の内部に設けられた処理部が位置算出処理を行うものとして説明した。つまり、上記の実施例では、ＧＰＳ受信部（ベースバンド処理回路部）が位置算出装置として機能するものとして説明した。しかし、位置算出の主体を、電子機器のプロセッサーであるホスト処理部としてもよい。
【０２３４】
この場合は、例えば、ベースバンド回路部の処理部がＧＰＳ衛星信号を捕捉してメジャメント情報を算出・取得し、ホスト処理部に出力する。そして、ホスト処理部は、処理部から入力したメジャメント情報を用いて、上記の原理に基づく位置算出演算を行う。この場合は、電子機器（ホスト処理部）が位置算出装置として機能することになる。
【符号の説明】
【０２３５】
１ 携帯型電話機、 ５ ＧＰＳアンテナ、 １０ ＧＰＳ受信部、 １１ ＲＦ受信回路部、 ２０ ベースバンド処理回路部、 ２１ 処理部、 ２３ 記憶部、 ３０ ホスト処理部、 ４０ 操作部、 ５０ 表示部、 ６０ 携帯電話用アンテナ、 ７０ 携帯電話用無線通信回路部、 ８０ 記憶部、 ９０ 時計部
【特許請求の範囲】
【請求項１】
測位用信号を受信することと、
前記測位用信号の受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行うことと、
を含む位置算出方法。
【請求項２】
過去の前記位置算出演算の結果に基づいて前記変動の基準を更新することと、
前記基準を用いて前記変動を算出することと、
を更に含む、
請求項１に記載の位置算出方法。
【請求項３】
前記位置算出演算は、少なくとも前記演算結果の変動を確率変数とした正規分布モデルに基づいて定められた演算であり、
過去の前記位置算出演算の結果に基づいて前記正規分布モデルを更新することを更に含む、
請求項１又は２に記載の位置算出方法。
【請求項４】
前記受信信号に基づいて擬似距離を測定することを更に含み、
前記位置算出演算は、前記演算結果の変動を確率変数とした正規分布モデルと、前記擬似距離に含まれる誤差を確率変数とした確率分布モデルとに基づいて定められた演算である、
請求項３に記載の位置算出方法。
【請求項５】
前記擬似距離に含まれる誤差を確率変数とした確率分布モデルは、正規分布モデル又は正規混合分布モデルである、
請求項４に記載の位置算出方法。
【請求項６】
前記位置算出演算は、位置情報の他に、クロックバイアス、速度情報及びクロックドリフトのうちの何れかを算出する演算であり、
前記正規分布モデルは、当該演算結果の各要素の変動を確率変数としたモデルである、
請求項３〜５の何れか一項に記載の位置算出方法。
【請求項７】
前記位置算出演算によって算出可能な算出パラメーターのうちの所定のパラメーターの値を、前記受信信号を用いた所定の収束演算を行って算出することを更に含み、
前記位置算出演算は、過去に前記収束演算によって算出された前記所定のパラメーターの値を用いる演算であることを含む、
請求項１〜６の何れか一項に記載の位置算出方法。
【請求項８】
測位用信号を受信する受信部と、
前記受信部による受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う位置算出部と、
を備えた位置算出装置。
【請求項１】
測位用信号を受信することと、
前記測位用信号の受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行うことと、
を含む位置算出方法。
【請求項２】
過去の前記位置算出演算の結果に基づいて前記変動の基準を更新することと、
前記基準を用いて前記変動を算出することと、
を更に含む、
請求項１に記載の位置算出方法。
【請求項３】
前記位置算出演算は、少なくとも前記演算結果の変動を確率変数とした正規分布モデルに基づいて定められた演算であり、
過去の前記位置算出演算の結果に基づいて前記正規分布モデルを更新することを更に含む、
請求項１又は２に記載の位置算出方法。
【請求項４】
前記受信信号に基づいて擬似距離を測定することを更に含み、
前記位置算出演算は、前記演算結果の変動を確率変数とした正規分布モデルと、前記擬似距離に含まれる誤差を確率変数とした確率分布モデルとに基づいて定められた演算である、
請求項３に記載の位置算出方法。
【請求項５】
前記擬似距離に含まれる誤差を確率変数とした確率分布モデルは、正規分布モデル又は正規混合分布モデルである、
請求項４に記載の位置算出方法。
【請求項６】
前記位置算出演算は、位置情報の他に、クロックバイアス、速度情報及びクロックドリフトのうちの何れかを算出する演算であり、
前記正規分布モデルは、当該演算結果の各要素の変動を確率変数としたモデルである、
請求項３〜５の何れか一項に記載の位置算出方法。
【請求項７】
前記位置算出演算によって算出可能な算出パラメーターのうちの所定のパラメーターの値を、前記受信信号を用いた所定の収束演算を行って算出することを更に含み、
前記位置算出演算は、過去に前記収束演算によって算出された前記所定のパラメーターの値を用いる演算であることを含む、
請求項１〜６の何れか一項に記載の位置算出方法。
【請求項８】
測位用信号を受信する受信部と、
前記受信部による受信信号を用いた演算であって、少なくとも演算結果の変動を確率変数とした所与の確率分布モデルに基づいて定められた位置算出演算を行う位置算出部と、
を備えた位置算出装置。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【公開番号】特開２０１２−２２０３７１（Ｐ２０１２−２２０３７１Ａ）
【公開日】平成２４年１１月１２日（２０１２．１１．１２）
【国際特許分類】
【出願番号】特願２０１１−８７４２２（Ｐ２０１１−８７４２２）
【出願日】平成２３年４月１１日（２０１１．４．１１）
【出願人】（０００００２３６９）セイコーエプソン株式会社 (51,324)
【Ｆターム（参考）】
【公開日】平成２４年１１月１２日（２０１２．１１．１２）
【国際特許分類】
【出願日】平成２３年４月１１日（２０１１．４．１１）
【出願人】（０００００２３６９）セイコーエプソン株式会社 (51,324)
【Ｆターム（参考）】
[ Back to top ]