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