(57)【特許請求の範囲】
【請求項１】
複数の測位衛星から送信される測位用信号を用いて得られる前記複数の測位衛星のそれぞれと受信機との距離から該受信機の位置を測位する単独測位装置において、
前記測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値から前記測位衛星の軌道情報および衛星時計誤差を観測する衛星情報検出手段と、
電離層遅延量情報を取得する電離層遅延情報取得手段と、
対流圏遅延量情報を取得する対流圏遅延情報取得手段と、
受信機位置を過去の受信機位置推定結果および前記測位衛星の軌道情報を用いて線形近似化し、
未知数である線形近似された受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とする回帰方程式を構成し、該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算する測位演算手段と、
を備えたことを特徴とする単独測位装置。
【請求項２】
前記測位演算手段は、整数値バイアス推定法を用いて前記整数値バイアスをフィックスして、前記受信機位置を推定演算する請求項１に記載の単独測位装置。
【請求項３】
複数エポックで取得した前記キャリア位相、前記コード擬似距離、前記受信機時計誤差、前記衛星軌道情報、前記衛星時計誤差、前記電離層遅延量情報、および、前記対流圏遅延量情報を記憶する一時記憶手段を備え、
前記測位演算手段は、取得した前記複数エポックの各情報を用いて、前記回帰方程式を構成する請求項１または請求項２に記載の単独測位装置。
【請求項４】
前記測位演算手段は、前記回帰方程式に状態推定アルゴリズムを適用する請求項１〜３のいずれかに記載の単独測位装置。
【請求項５】
前記測位演算手段は、前記衛星軌道情報を目的変数に含み、前記衛星軌道情報に対応する衛星軌道誤差を説明変数に含む前記回帰方程式を用いて、測位衛星位置を推定演算する請求項１〜４のいずれかに記載の単独測位装置。
【請求項６】
複数の測位衛星から送信される測位用信号を用いて得られる前記複数の測位衛星のそれぞれと受信機との距離から該受信機の位置を測位する単独測位方法であって、
前記測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値から前記測位衛星の軌道情報および衛星時計誤差を観測し、
電離層遅延量情報、および対流圏遅延量情報を取得する対流圏遅延情報取得し、
過去の受信機位置推定結果および前記測位衛星の軌道情報を用いて受信機位置を線形近似化し、
未知数である線形近似された受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とする回帰方程式を構成し、
該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算する、
ことを特徴とする単独測位方法。
【請求項７】
整数値バイアス推定法を用いて前記整数値バイアスをフィックスして、前記受信機位置を推定演算する請求項６に記載の単独測位方法。
【請求項８】
複数エポックで取得した前記キャリア位相、前記コード擬似距離、前記受信機時計誤差、前記衛星軌道情報、前記衛星時計誤差、前記電離層遅延量情報、および、前記対流圏遅延量情報を一時記憶し、
取得した前記複数エポックの各情報を用いて、前記回帰方程式を構成する請求項６または請求項７に記載の単独測位方法。
【請求項９】
前記回帰方程式に状態推定アルゴリズムを適用する請求項６〜８のいずれかに記載の単独測位方法。
【請求項１０】
前記衛星軌道情報を目的変数に含み、前記衛星軌道情報に対応する衛星軌道誤差を説明変数に含む前記回帰方程式を用いて、測位衛星位置を推定演算する請求項６〜９のいずれかに記載の単独測位方法。
【請求項１１】
複数の測位衛星から送信される測位用信号を用いて得られる前記複数の測位衛星のそれぞれと受信機との距離から該受信機の位置を測位する単独測位装置において、
前記測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値から前記測位衛星の軌道情報および衛星時計誤差を観測する衛星情報検出手段と、
電離層遅延量情報を取得する電離層遅延情報取得手段と、
対流圏遅延量情報を取得する対流圏遅延情報取得手段と、
過去の受信機位置推定結果および前記測位衛星の軌道情報を用いて、前記受信機と前記測位衛星との距離を推定受信機位置の周りで一次ティラー級数展開により線形近似することで、前記受信機位置を線形近似化して設定し、
未知数である受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、
観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とし、
前記一次ティラー級数展開で得られる前記距離の線形近似値を受信機位置の各座標で偏微分してなる係数を演算行列要素として用いることで、回帰方程式を構成し、
該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算する測位演算手段と、
を備えたことを特徴とする単独測位装置。
【請求項１２】
複数の測位衛星から送信される測位用信号を用いて得られる前記複数の測位衛星のそれぞれと受信機との距離から該受信機の位置を測位する単独測位方法であって、
前記測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値から前記測位衛星の軌道情報および衛星時計誤差を観測し、
電離層遅延量情報、および対流圏遅延量情報を取得する対流圏遅延情報取得し、
過去の受信機位置推定結果および前記測位衛星の軌道情報を用いて、前記受信機と前記測位衛星との距離を推定受信機位置の周りで一次ティラー級数展開により線形近似することで、前記受信機位置を線形近似化して設定し、
未知数である受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、
観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とし、
前記一次ティラー級数展開で得られる前記距離の線形近似値を受信機位置の各座標で偏微分してなる係数を演算行列要素として用いることで、回帰方程式を構成し、
該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算する、
ことを特徴とする単独測位方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、測位用信号を用いて受信機の単独測位を行う単独測位装置および単独測位方法に関するものである。
【背景技術】
【０００２】
従来、測位衛星から送信される測位用信号を用いて単独測位を行う装置が各種開示されている。これらの装置による基本測位演算は、コード擬似距離を用いて、受信機の３次元の位置（誤差）と受信機の時計誤差とを未知数とした非線形連立方程式にニュートン法や拡張カルマンフィルタを適用することで求めるものであった。この際、電離層遅延による誤差や対流圏遅延による誤差の影響を除去する方法が各種考案されている。その１つとして、これらの遅延による影響を初期状態から「０」に設定する方法があり、またその１つとして、１重位相差等の演算処理を行ってこれらの遅延量を推定する方法があった。
【非特許文献１】日本測地学会編著，「新訂版ＧＰＳ−人工衛星による精密測位システム−」，社団法人日本測量協会，１９８９年１１月１５日，ｐ．１２１−１４０
【発明の開示】
【発明が解決しようとする課題】
【０００３】
ところが、従来のコード擬似距離のみを用いた単独測位装置および単独測位方法では、現実に存在する電離層遅延や対流圏遅延を無視したり、あるいは十分な精度で求めずして推定演算を行っていた。このため、測位の推定値のばらつきが大きく、高精度な測位を行うことができなかった。また、コード擬似距離を用いるとともに、電離層遅延情報や対流圏遅延情報を用いて１重位相差処理等を行う場合、前述のコード擬似距離のみを用いる場合よりも推定精度は向上するものの、推定演算内の１重位相差処理を行う等、推定演算処理が複雑になったり、１重位相差処理を行った後のノイズの設定により大きく推定演算結果が異なったりする。このため、推定演算の処理内容が複雑なわりには、推定演算結果の精度の向上性が悪く、さらに推定演算処理速度が遅くなる可能性があった。
【０００４】
したがって、この発明の目的は、推定演算処理を複雑にすることなく確実に高精度測位を実現することができる単独測位装置および単独測位方法を提供することにある。
【課題を解決するための手段】
【０００５】
この発明は、複数の測位衛星から送信される測位用信号を用いて得られる複数の測位衛星のそれぞれと受信機との距離から該受信機の位置を測位する単独測位装置において、
測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値から測位衛星の軌道情報および衛星時計誤差を観測する衛星情報検出手段と、電離層遅延量情報を取得する電離層遅延情報取得手段と、対流圏遅延量情報を取得する対流圏遅延情報取得手段と、受信機位置を過去の受信機位置推定結果および測位衛星の軌道情報を用いて線形近似化し、未知数である線形近似された受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とする回帰方程式を構成し、該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算する測位演算手段と、を備えたことを特徴としている。
【０００６】
また、この発明は、複数の測位衛星から送信される測位用信号を用いて得られる複数の測位衛星のそれぞれと受信機との距離から該受信機の位置を測位する単独測位方法であって、
測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値から測位衛星の軌道情報および衛星時計誤差を観測し、電離層遅延量情報、および対流圏遅延量情報を取得する対流圏遅延情報取得し、過去の受信機位置推定結果および測位衛星の軌道情報を用いて受信機位置を線形近似化し、未知数である線形近似された受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とする回帰方程式を構成し、
該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算することを特徴としている。
【０００７】
これらの構成では、搬送波毎で測位衛星毎のキャリア位相、ＰＮコード毎で測位衛星毎のコード擬似距離、測位衛星毎の衛星軌道情報、時計誤差、電離層遅延情報、対流圏遅延情報、を観測値とし、受信機の３次元位置、受信機の時計誤差、搬送波毎で測位衛星毎の整数値バイアスを未知数する。そして、各観測値を用いて目的変数を形成し、各未知数を用いて説明変数を形成することで回帰方程式が構成される。この際、受信機位置は、過去の受信機位置推定結果および測位衛星の軌道情報を用いて線形近似される。そして、この回帰方程式に最小二乗法等のパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を含む未知数が推定演算される。
【０００８】
また、この発明の単独測位装置および単独測位方法は、整数値バイアス推定法を用いて整数値バイアスをフィックスして受信機位置を推定演算することを特徴としている。
【０００９】
この構成では、搬送波の整数値バイアスに対してＬＡＭＢＤＡ法等の整数値バイアス推定法を適用することで、整数値バイアスが整数値としてフィックスされる。このフィックスされた整数値バイアスを回帰方程式に既知数として代入することで、回帰方程式内の未知数の数が減少し、推定精度が向上する。
【００１０】
また、この発明の単独測位装置および単独測位方法は、複数エポックで取得したキャリア位相、コード擬似距離、受信機時計誤差、衛星軌道情報、衛星時計誤差、電離層遅延量情報、および、対流圏遅延量情報を一時記憶手段で記憶し、測位演算手段で取得した複数エポックの各情報を用いて回帰方程式を構成することを特徴としている。
【００１１】
この構成では、複数エポックに亘る各情報を得ることで、観測値の数と未知数の数とが増加するが、観測値の数の増加の方が多いので、回帰方程式の未知数を推定がより容易になる。
【００１２】
また、この発明の単独測位装置および単独測位方法は、測位演算手段で、回帰方程式に状態推定アルゴリズムを適用することを特徴としている。
【００１３】
この構成では、前述の回帰方程式にカルマンフィルタや非線形フィルタ等の状態推定アルゴリズムを適用することで、移動する受信機位置の推定が可能となる。
【００１４】
また、この発明の単独測位装置および単独測位方法は、測位演算手段で、衛星軌道情報を目的変数に含み、衛星軌道誤差を説明変数に含む回帰方程式を用いて、測位衛星位置を推定演算することを特徴としている。
【００１５】
この構成では、測位衛星位置に対応する衛星軌道誤差を未知数として説明変数に加えることで、回帰方程式から測位衛星位置が推定演算される。
【発明の効果】
【００１６】
この発明によれば、電離層遅延量情報、対流圏遅延量情報、衛星軌道情報、衛星時計誤差を観測値として用い、且つキャリア位相およびコード擬似距離も観測値として用いて、受信機位置、整数値バイアスを未知数とする回帰方程式に適用することで、受信機位置および整数値バイアスの推定精度が向上する。これにより、１つの回帰方程式のみを用いるという簡素な演算で、受信機の単独測位を高精度に行うことができる。
【００１７】
また、この発明によれば、ＬＡＭＢＤＡ法等の整数値バイアス推定法を用いることで、より高精度に単独測位を行うことができる。
【００１８】
また、この発明によれば、複数エポックに亘る観測値を用いることで必要な観測衛星数を少なくして、高精度な単独測位を行うことができる。
【００１９】
また、この発明によれば、前述の線形回帰方程式にカルマンフィルタ等の状態推定アルゴリズムを適用することで、移動する受信機の位置を高精度に測位することができる。
【００２０】
また、この発明によれば、衛星軌道情報を観測値とし、衛星軌道誤差を未知数として推定演算することで、受信機と測位衛星との間の距離をより高精度に測位することができる。これにより、受信機位置をより高精度に測位することができる。
【発明を実施するための最良の形態】
【００２１】
本発明の実施形態に係る単独測位装置について図を参照して説明する。なお、以下の説明では、ＧＰＳについて説明するが、他の全てのＧＮＳＳ（全地球的航法衛星システム）に適用することができる。
図１は本実施形態の単独測位装置の概略構成を示すブロック図である。
また、図２は本実施形態の単独測位装置とＧＰＳ受信機とからなる単独測位システムの処理フローを示すフローチャートである。
【００２２】
図１に示すように、単独測位装置１０は、ＧＰＳ受信機３０に接続し、航法メッセージ解析部１１、衛星情報処理部１２、および測位演算部１３を備える。
【００２３】
ＧＰＳ受信機３０はアンテナ２０に接続し、アンテナ２０で受信したＧＰＳ衛星（測位衛星）からのＧＰＳ信号（測位用信号）より、既知の方法でＬ１波およびＬ２波のキャリア位相、Ｃ／ＡコードおよびＰ（Ｙ）コードのコード擬似距離（擬似距離）を取得するともに、Ｌ１波に重畳された航法メッセージを取得する（ｓ１→ｓ２）。そして、ＧＰＳ受信機３０は航法メッセージを航法メッセージ解析部１１に出力する。また、ＧＰＳ受信機３０はキャリア位相およびコード擬似距離を測位演算部１３に出力する。
【００２４】
航法メッセージ解析部１１は、入力された航法メッセージを解析して、電離層遅延量情報、各ＧＰＳ衛星の時計誤差、軌道情報を取得する。また、方法メッセージ解析部１１は、数式モデルを用いて対流圏遅延量情報を取得する。そして、航法メッセージ解析部１１は取得した各情報を衛星情報処理部１２に出力する。
【００２５】
衛星情報処理部１２は、ＧＰＳ衛星のエフェメリス情報を用いて測位に用いるＧＰＳ衛星を選定して、選定したＧＰＳ衛星に関する衛星軌道情報、衛星時計誤差情報、電離層遅延量情報、および対流圏量遅延情報を測位演算部１３に出力する。
【００２６】
測位演算部１３は、入力された各観測値、すなわち、選定された各ＧＰＳ衛星に対するＬ１波のキャリア位相、Ｌ２波のキャリア位相、Ｃ／Ａコードによるコード擬似距離、Ｐ（Ｙ）コードのコード擬似距離、各ＧＰＳ衛星の衛星軌道情報、各ＧＰＳ衛星の衛星時計誤差情報、電離層遅延量情報、対流圏遅延量情報を用いて、後述する線形回帰方程式を形成する。そして、測位演算部１３は、この線形回帰方程式にパラメータ推定アルゴリズムである最小二乗法を適用させて、Ｌ１、Ｌ２波に対する整数値バイアスＮL1，ＮL2と受信機の位置とを推定演算する（ｓ４）。この演算は推定値の変化が予め設定された所定の閾値以下になるまで繰り返し行われ、推定値の変化が所定閾値に達した時点での受信機位置の推定演算結果が測位結果として出力される。
【００２７】
ここで、航法メッセージ解析部１１、衛星情報処理部１２、測位演算部１３は、それぞれ、以下に示すアルゴリズムを実現するマイクロプロセッサ等の数値演算処理器からなる。そして、これらの部分は複数の数値演算処理器により形成してもよく、１つの数値演算処理器により形成してもよい。
【００２８】
次に、前述の整数値バイアスＮL1，ＮL2および受信機位置ｕの推定演算アルゴリズムについて詳述する。
【００２９】
一般に、受信機ｕ、ＧＰＳ衛星ｐに対するキャリア位相φpL,uの観測方程式は式（１）で表すことができ、コード擬似距離（擬似距離）ρpc,uの観測方程式は式（２）で表される。ここで、マルチパス誤差は微少として無視する。
【００３０】
【数１】
0004146877.tif 000002
【００３１】
ここで、λLはＬ波の波長を示し、ｒpu（ｔ，ｔ−τpu）は時刻ｔでの受信機ｕと時刻（ｔ−τpu）でのＧＰＳ衛星との距離を示す。δＩpu（ｔ）はＬ１波の電離層遅延量を示し、δＴpu（ｔ）はＬ１波、Ｌ２波の対流圏量遅延を示す。δｔu（ｔ）は真の時刻ｔでの受信機ｕの時計誤差を示し、δｔp（ｔ−τpu）は時刻（ｔ−τpu）でのＧＰＳ衛星ｐの時計誤差を示す。ＮpL,uはＬ波における受信機ｕとＧＰＳ衛星ｐとの間の整数値バイアスを示し、εpL,u（ｔ），ｅpc,u（ｔ）はそれぞれ観測雑音を示す。
【００３２】
したがって、Ｌ１波におけるキャリア位相φpL1,uの観測方程式は式（３）となり、Ｌ２波におけるキャリア位相φpL2,uの観測方程式は式（４）となる。
【００３３】
【数２】
0004146877.tif 000003
【００３４】
ここで、ｆL1はＬ１波の周波数、ｆL2はＬ２の周波数を示す。
【００３５】
さらに、ΦpL1,u，ΦpL2,uをそれぞれ次式（３’），（４’）により定義する。
【００３６】
【数３】
0004146877.tif 000004
【００３７】
また、Ｃ／Ａコードによるコード擬似距離ρpCA,uの観測方程式は式（５）となり、Ｐコードによるコード擬似距離ρpP,uの観測方程式は式（６）となる。
【００３８】
【数４】
0004146877.tif 000005
【００３９】
ところで、受信機とＧＰＳ衛星との距離ｒpu（ｔ，ｔ−τpu）は、式（７）で表すことができる。
【００４０】
【数５】
0004146877.tif 000006
【００４１】
次に、未知数である受信機位置ｕ（ｔ）≡［ｘu（ｔ），ｙu（ｔ），ｚu（ｔ）］Tを先験的な推定受信機位置ｕ(j)（ｔ）≡［ｘu(j)（ｔ），ｙu(j)（ｔ），ｚu(j)（ｔ）］Tのまわりで１次のティラー級数展開を行い、ｒpu（ｔ）を線形近似すると次式を得られる。
【００４２】
【数６】
0004146877.tif 000007
【００４３】
この式（９）により、式（３’）、（４’）、（５）、（６）は、次式（１０）〜（１３）で表される。
【００４４】
【数７】
0004146877.tif 000008
【００４５】
ここで、Φep,(j)L1,u，Φep,(j)L2,u，ρep,(j)L1,u，ρep,(j)L2,u、を式（１４）〜（１７）で定義する。
【００４６】
【数８】
0004146877.tif 000009
【００４７】
これにより、式（１０）〜（１３）は次式（１８）〜（２１）となる。
【００４８】
【数９】
0004146877.tif 000010
【００４９】
これは、すなわち、キャリア位相、コード擬似距離を目的変数とし、受信機の位置、電離層遅延量、対流圏遅延量、整数値バイアス、誤差要因を説明変数とする近似的な線形回帰方程式に相当する。
【００５０】
ここで、Ｇ(j)を式（２２）で定義する。
【００５１】
【数１０】
0004146877.tif 000011
【００５２】
なお、Ｇ(j)uは式（２３）で表される。
【００５３】
【数１１】
0004146877.tif 000012
【００５４】
また、Φep,(j)L1,u，Φep,(j)L2,u，ρep,(j)L1,u，ρep,(j)L2,u、を行列式で表す。
【００５５】
【数１２】
0004146877.tif 000013
【００５６】
この式を用いることで、式（１８）〜式（２１）はベクトル行列からなる次式で表される。
【００５７】
【数１３】
0004146877.tif 000014
【００５８】
もしくは、
【００５９】
【数１４】
0004146877.tif 000015
【００６０】
ところで、前述のように、航法メッセージにはＧＰＳ衛星毎の時計誤差情報が含まれており、航法メッセージ解析部１１でＧＰＳ衛星毎の時計誤差情報が取得されている。この衛星時計誤差の観測値δｔｅpは次式で表され、前述の回帰方程式に加えることができる。
【００６１】
【数１５】
0004146877.tif 000016
【００６２】
また、航法メッセージから得られるＧＰＳ衛星の情報からＧＰＳ衛星毎の電離層遅延量および対流圏遅延量をそれぞれ電離層遅延モデルおよび対流圏遅延モデルを用いて取得することができる。これにより、電離層遅延量の観測値δＩｅuとし、対流圏遅延量の観測値δＴｅuはそれぞれ次式で表され、前述の回帰方程式に加えることができる。
【００６３】
【数１６】
0004146877.tif 000017
【００６４】
この結果、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量、対流圏遅延量を目的変数とし、整数値バイアス、受信機位置を説明変数として含む線形回帰方程式が実現される。
【００６５】
すなわち、式（２６）に示した回帰方程式は、次式に置き換えられる。
【００６６】
【数１７】
0004146877.tif 000018
【００６７】
ここで、υuEの共分散行列をＲは式（３５）で表される。
【００６８】
【数１８】
0004146877.tif 000019
【００６９】
従って、θ(j+1)の推定値θｅ(j+1)は、式（３６）となる。
【００７０】
【数１９】
0004146877.tif 000020
【００７１】
また、推定値θｅ(j+1)の分散値は、式（３７）となる。
【００７２】
【数２０】
0004146877.tif 000021
【００７３】
この式（３６）に基づいて、最小二乗法を用いて繰り返し推定演算を行うことで、θを構成する未知数を収束させて各値を取得する。この際、説明変数ベクトルθを構成する未知数の収束条件は、未知数の推定値の繰り返し演算での変化量が、予め設定した閾値以下に達した時点である。そして、この条件下で説明変数ベクトルθの構成要素の値を確定させる。このような演算を行うことで、整数値バイアスＮL1，ＮL2や受信機位置ｕを推定することができる。
【００７４】
このように、本実施形態の構成を用いることで、１重位相差等の処理を行うことなく整数値バイアスＮL1，ＮL2や受信機位置ｕの推定演算、すなわち単独測位を行うことができる。これにより、従来例より簡素なアルゴリズムで単独測位を行うことができ、且つ単独測位の演算処理速度を向上することが可能になる。
【００７５】
ところで、このような推定演算を用いる場合、未知パラメータはθuの各要素であり、式（３１）に示す回帰方程式内の未知パラメータ数は、受信機位置のｕが３個、受信機時計誤差δｔuが１個、ＧＰＳ衛星の時計誤差δｔpがｎs個（ＧＰＳ衛星数）、電離層遅延δＩがｎs個、対流圏遅延δＴがｎs個、Ｌ１波整数値バイアスＮL1がｎs個、Ｌ２波整数値バイアスＮL2がｎs個、の合計４＋５ｎs個となる。一方で、観測値は、Ｌ１波のキャリア位相Φｅ(j)L1,u、Ｌ２波のキャリア位相Φｅ(j)L2,u、ＣＡコード擬似距離ρｅ(j)CA,u、Ｐコード擬似距離ρｅ(j)CA,u、ＧＰＳ衛星の時計誤差δｔｅp、電離層遅延δＩｅ、対流圏遅延δＴｅがそれぞれ測位可能なＧＰＳ衛星数に相当するｎsであり、合計７ｎs個となる。
【００７６】
従って、この回帰方程式で前記未知数の解を確定する成立条件は、次式となる。
【００７７】
４＋５ｎs≦７ｎs
これはｎs≧２を意味する。すなわち、少なくとも２つのＧＰＳ衛星からの測位用信号を受信すれば受信機の単独測位が可能となる。しかしながら、この回帰方程式を最小二乗法で解く場合、式（３３）でのＨ(j)Eの逆行列が必要となる。ところが、測位衛星数ｎsが２の場合には前記行列が特異行列となり、ｎsが３の場合には前記行列が特異となる可能性が非常に高くなり、逆行列が存在せず現実的に演算できない。
【００７８】
また、前述の実施形態では、Ｐコードを用いて推定演算を行ったが、Ｐ（Ｙ）コードは秘匿コードであるため、コード擬似距離の観測は事実上難しい。この場合、観測値が１種類（ｎs個）減少するので回帰方程式の成立条件は次式となる。
【００７９】
４＋５ｎs≦６ｎs
これはｎs≧４を意味し、実質的に少なくとも４つの測位衛星からの測位用信号を受信すれば受信機の単独測位が可能となる。すなわち、少なくとも４つの測位衛星からの測位用信号を受信できれば、前述の単独測位は可能となる。
【００８０】
次に、本実施形態の単独測位方法を用いた場合のシミュレーション結果について説明する。なお、このシミュレーション結果は、Ｐコードを用いずに推定演算を行ったものである。
【００８１】
図３は本実施形態の単独測位方法を用いた場合の受信機位置（経度、緯度）を示した図である。図４は本実施形態の単独測位方法を用いた場合の受信機位置（楕円体高）を示した図である。
【００８２】
図３、図４において、円形プロットは本実施形態の単独測位結果１０１を示す。四角形プロットは従来のコード擬似距離のみを用いた単独測位結果１０２を示す。三角形プロットは従来のコード擬似距離、電離層遅延量、対流圏遅延量、および衛星時計誤差を用いた単独測位結果１０３を示す。また、図３におけるアスタリスクマーキング点および図４における実線は相対測位結果を示す。
【００８３】
また、表１は本実施形態の単独測位方法による受信機位置（経度、緯度、楕円体高）の平均および標準偏差を示す。表２は、従来のコード擬似距離のみを用いた単独測位方法による受信機位置（経度、緯度、楕円体高）の平均および標準偏差を示す。表３は、従来のコード擬似距離、電離層遅延量、対流圏遅延量、衛星時計誤差を用いた単独測位方法による受信機位置（経度、緯度、楕円体高）の平均および標準偏差を示す。
【００８４】
なお、このシミュレーションでは、受信機位置の初期値として、ＲＩＮＥＸデータの”ＡＰＰＲＯＸ  ＰＯＳＩＴＩＯＮ  ＸＹＺ”に記載された座標を用いた。また、電離層遅延量は、航法メッセージに含まれる電離層遅延に関するデータから、いわゆる放送モデル（Ｋｌｏｂｕｃｈａｒモデル）を適用して算出した値を用い、対流圏遅延量は、航法データから計算される測位衛星の仰角情報から次式を適用して算出した値を用いている。
【００８５】
【数２１】
0004146877.tif 000022
【００８６】
ここで、ξは仰角である。
さらに、観測誤差の分散は、それぞれ、コード擬似距離が１．５［ｍ］（１σ）で、キャリア位相が（波長λ／１０）＋（１．５／１０）［ｍ］（１σ）で、衛星時計誤差が３．６［ｍ］（１σ）で、電離層遅延が７．０［ｍ］（１σ）で、対流圏遅延が０．７［ｍ］（１σ）でそれぞれ設定されている。また、このシミュレーションは、最小二乗法による繰り返し演算での、受信機位置の推定値の変化のノルム値が１×１０-3［ｍ］以下になるまで計算を繰り返した。
【００８７】
【表１】
0004146877.tif 000023
【００８８】
【表２】
0004146877.tif 000024
【００８９】
【表３】
0004146877.tif 000025
【００９０】
これらの結果に示すように、本実施形態の単独測位方法を用いることにより、高精度に受信機位置を測位することができる。
【００９１】
以上のように、本実施形態の単独測位装置および単独測位方法を用いることにより、従来よりも簡素な演算処理で高精度に受信機の単独測位を行うことができる。
【００９２】
なお、前述の最小二乗法による推定演算時に、式（３７）より得られる分散値を用いて、ＬＡＭＢＤＡ法を適用することで、整数値バイアスＮL1，ＮL2をフィックスして、さらに高精度の受信機位置を測位することができる。
【００９３】
また、前述の説明では、１エポックの観測値を用いて単独測位を行う方法を示したが複数エポックの観測値を用いても単独測位を行うことができる。この際、単独測位装置には、複数エポックのキャリア位相、コード擬似距離、測位衛星の衛星軌道情報、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を記憶するメモリ（本発明の「一時記憶手段」に相当する。）を備える。
【００９４】
このような場合、例えば、Ｃ／Ａコードのみを使用した場合、Ｋエポックで取得できる観測値は、次式で表される。
【００９５】
６ｎs×Ｋ＝６Ｋｎs
一方、観測エポックが複数となるキネマティック測位では、整数値バイアスＮL1，ＮL2を除く各未知数（未知パラメータ）がエポック数に応じて増加するので、未知数の数は、次式で表される。
【００９６】
４＋５ｎs＋（Ｋ−１）（４＋３ｎs）＝４Ｋ＋（３Ｋ＋２）ｎs
したがって、回帰方程式を最小二乗法により解くことができる必要条件は、次式となる。
【００９７】
６Ｋｎs≧４Ｋ＋（３Ｋ＋２）ｎs
（３Ｋ−２）ｎs≧４Ｋ
すなわち、
ｎs≧４Ｋ／（３Ｋ−２）
そして、エポック数Ｋを２→∞の場合、ｎsの取り得る整数値を計算すると、ＧＰＳ衛星数ｎsは次式の値となる。
【００９８】
ｎs≧２
これにより、観測エポック数を増加させることにより、観測するＧＰＳ衛星数を最低２個として測位を行なうことができる。
【００９９】
ここで、このようなキネマティック測位において、カルマンフィルタを適用することで、より高精度な測位が可能となる。以下に、カルマンフィルタを適用した場合の測位アルゴリズムについて説明する。
【０１００】
まず、式（３１）より、エポック数に対応する観測時間ｔを要素として表現し直すと、
【０１０１】
【数２２】
0004146877.tif 000026
【０１０２】
これを、ベクトル行列表示すると、式（３９）になる。
【０１０３】
【数２３】
0004146877.tif 000027
【０１０４】
ここでｃδｔp，δＩu，δＴuをそれぞれ式（４０ａ），（４０ｂ），（４０ｃ）とする。
【０１０５】
【数２４】
0004146877.tif 000028
【０１０６】
これらを式（３９）の回帰方程式に代入し、さらに、次式に示す定義を行う。
【０１０７】
【数２５】
0004146877.tif 000029
【０１０８】
この結果、新たな回帰方程式は、式（４５）で表される。
【０１０９】
【数２６】
0004146877.tif 000030
【０１１０】
ここで、未知の受信機位置ｕ（ｔ）の速度に対するマルコフ過程モデルを仮定し、かつ、受信機の時計誤差ｃδｔu（ｔ）に対してマルコフ過程モデルを仮定できるものとし、速度ベクトルをｖ（ｔ）とおき、受信機の新たな状態ベクトルηu（ｔ）を次式で定義する。
【０１１１】
【数２７】
0004146877.tif 000031
【０１１２】
そして、この状態ベクトルηu（ｔ）について次の状態方程式を構成する。
【０１１３】
【数２８】
0004146877.tif 000032
【０１１４】
この際、式（４５）の観測方程式は、次式で表すことができる。
【０１１５】
【数２９】
0004146877.tif 000033
【０１１６】
この観測方程式をベクトル表示すると、式（４９）となる。
【０１１７】
【数３０】
0004146877.tif 000034
【０１１８】
これは、すなわち、式（４７）、式（４８）（式（４９））からなるカルマンフィルタを意味する。この際、観測誤差υu,R（ｔ）の共分散行列Ｒは、式（３５）と同様に求めることができる。そして、これら状態方程式（４７）、観測方程式（４８）、を用いることにより移動する受信機位置を推定演算することができる。
【０１１９】
このような構成とすることで、受信機位置が移動してもカルマンフィルタにより受信機の移動を推定しながら、受信機位置を高精度に測位することができる。
【０１２０】
ところで、前述の説明では、ＧＰＳ衛星の位置を推定演算しなかったが、次に示す方法（アルゴリズム）を用いることにより、ＧＰＳ衛星の位置を推定演算することができる。
【０１２１】
式（７）に示すように、受信機とＧＰＳ衛星との距離を定義し、受信機とともにＧＰＳ衛星の位置を推定演算する場合、受信機位置に関する線形近似式は、前述のように表される。
【０１２２】
【数３１】
0004146877.tif 000035
【０１２３】
測位衛星位置に関する線形近似式は、次式で表される。
【０１２４】
【数３２】
0004146877.tif 000036
【０１２５】
そして、ｓｐ≡［ｘp，ｙp，ｚp］Tと定義すると次式の関係が成り立つ。
【０１２６】
【数３３】
0004146877.tif 000037
【０１２７】
したがって、ｕ，ｓｐの推定値は、それぞれｕ(j)，ｓｅpと１次ティラー級数展開とを用いて次式で近似される。
【０１２８】
【数３４】
0004146877.tif 000038
【０１２９】
したがって、式（１０）から式（１３）は、それぞれ次式で表される。
【０１３０】
【数３５】
0004146877.tif 000039
【０１３１】
ここで、Φep,(j)L1,u，Φep,(j)L2,u，ρep,(j)L1,u，ρep,(j)L2,u、をそれぞれ式（５５）〜（５８）で定義する。
【０１３２】
【数３６】
0004146877.tif 000040
【０１３３】
これにより、式（５１）〜（５４）は次式（５９）〜（６２）となる。
【０１３４】
【数３７】
0004146877.tif 000041
【０１３５】
ここで、Φep,(j)L1,u，Φep,(j)L2,u，ρep,(j)L1,u，ρep,(j)L2,u、をそれぞれ次式で定義する。
【０１３６】
【数３８】
0004146877.tif 000042
【０１３７】
これにより、式（５９）〜式（６２）はベクトル行列からなる次式で表される。
【０１３８】
【数３９】
0004146877.tif 000043
【０１３９】
また、精密軌道衛星に関する観測値ベクトルｓｋを次式で表す。
【０１４０】
【数４０】
0004146877.tif 000044
【０１４１】
これを式（６３）に適用すると式（６４）となる。
【０１４２】
【数４１】
0004146877.tif 000045
【０１４３】
このような回帰方程式を用いることで、ＧＰＳ衛星の位置についても推定演算することができる。そして、ＧＰＳ衛星位置の推定位置精度が向上することで、受信機位置の推定位置をより高精度に推定することができる。
【０１４４】
さらに、以下のアルゴズムを用いることで、このような場合でもカルマンフィルタを適用することができる。
【０１４５】
すなわち、ｓ，ｃδｔp，δＩu，δＴuをそれぞれ式（６５ａ）〜（６５ｄ）で定義する。
【０１４６】
【数４２】
0004146877.tif 000046
【０１４７】
これらの関係を用い、式（６４）の回帰方程式に代入し、さらに
次に示す定義を行う。
【０１４８】
【数４３】
0004146877.tif 000047
【０１４９】
この結果、新たな回帰方程式は、式（７０）で表される。
【０１５０】
【数４４】
0004146877.tif 000048
【０１５１】
そして、この観測方程式に対して前述のように状態方程式を定義することで、カルマンフィルタを構成することができる。
【０１５２】
なお、前述の実施形態では、線形回帰方程式に最小二乗法を適用した例を示したが、他のパラメータ推定アルゴリズムを用いてもよい。
【０１５３】
また、前述の実施形態では、受信機位置の線形近似に１次のティラー級数展開を用いたが、他の線形化演算を用いてもよい。
【０１５４】
また、前述の実施形態では、ＬＡＭＢＤＡ法を利用する例を示したが、他の整数値バイアス推定方法を用いてもよい。
【０１５５】
また、前述の実施形態では、カルマンフィルタを適用する例を示したが、他の状態推定アルゴリズムを用いてもよい。
【図面の簡単な説明】
【０１５６】
【図１】本発明の実施形態の単独測位装置の概略構成を示すブロック図である。
【図２】本発明の実施形態の単独測位装置とＧＰＳ受信機とからなる単独測位システムの処理フローを示すフローチャートである。
【図３】本発明の実施形態の単独測位方法を用いた場合の受信機位置（経度、緯度）を示した図である。
【図４】本発明の実施形態の単独測位方法を用いた場合の受信機位置（楕円体高）を示した図である。
【符号の説明】
【０１５７】
１０−単独測位装置
１１−航法メッセージ解析部
１２−衛星情報処理部
１３−測位演算部
２０−ＧＰＳアンテナ
３０−ＧＰＳ受信機
【図１】
0004146877.tif 000049
【図２】
0004146877.tif 000050
【図３】
0004146877.tif 000051
【図４】
0004146877.tif 000052