JPWO2006022318A1 - 単独測位装置および単独測位方法 - Google Patents

Abstract

入力された各観測値、すなわち、選定された各ＧＰＳ衛星に対するＬ１波のキャリア位相、Ｌ２波のキャリア位相、Ｃ／Ａコードのコード擬似距離、Ｐ（Ｙ）コードのコード擬似距離、時計誤差情報、電離層遅延情報、対流圏遅延情報を目的変数とし、整数値バイアス、および受信機位置を説明変数に含む線形回帰方程式を構成する。この際、受信機位置は、過去の受信機位置推定結果から線形近似しておく。そして、この線形回帰方程式に最小二乗法を適用させて、整数値バイアスと受信機位置とを推定演算する。

Description

本発明は、測位用信号を用いて受信機の単独測位を行う単独測位装置および単独測位方法に関するものである。

従来、測位衛星から送信される測位用信号を用いて単独測位を行う装置が各種開示されている。これらの装置による基本測位演算は、コード擬似距離を用いて、受信機の３次元の位置（誤差）と受信機の時計誤差とを未知数とした非線形連立方程式にニュートン法や拡張カルマンフィルタを適用することで求めるものであった。この際、電離層遅延による誤差や対流圏遅延による誤差の影響を除去する方法が各種考案されている。その１つとして、これらの遅延による影響を初期状態から「０」に設定する方法があり、またその１つとして、１重位相差等の演算処理を行ってこれらの遅延量を推定する方法があった。
日本測地学会編著，「新訂版ＧＰＳ−人工衛星による精密測位システム−」，社団法人日本測量協会，１９８９年１１月１５日，ｐ．１２１−１４０

ところが、従来のコード擬似距離のみを用いた単独測位装置および単独測位方法では、現実に存在する電離層遅延や対流圏遅延を無視したり、あるいは十分な精度で求めずして推定演算を行っていた。このため、測位の推定値のばらつきが大きく、高精度な測位を行うことができなかった。また、コード擬似距離を用いるとともに、電離層遅延情報や対流圏遅延情報を用いて１重位相差処理等を行う場合、前述のコード擬似距離のみを用いる場合よりも推定精度は向上するものの、推定演算内の１重位相差処理を行う等、推定演算処理が複雑になったり、１重位相差処理を行った後のノイズの設定により大きく推定演算結果が異なったりする。このため、推定演算の処理内容が複雑なわりには、推定演算結果の精度の向上性が悪く、さらに推定演算処理速度が遅くなる可能性があった。

したがって、この発明の目的は、推定演算処理を複雑にすることなく確実に高精度測位を実現することができる単独測位装置および単独測位方法を提供することにある。

この発明は、複数の測位衛星から送信される測位用信号を用いて得られる複数の測位衛星のそれぞれと受信機との距離から該受信機の位置を測位する単独測位装置において、
測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値から測位衛星の軌道情報および衛星時計誤差を観測する衛星情報検出手段と、電離層遅延量情報を取得する電離層遅延情報取得手段と、対流圏遅延量情報を取得する対流圏遅延情報取得手段と、受信機位置を過去の受信機位置推定結果および測位衛星の軌道情報を用いて線形近似化し、未知数である線形近似された受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とする回帰方程式を構成し、該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算する測位演算手段と、を備えたことを特徴としている。

また、この発明は、複数の測位衛星から送信される測位用信号を用いて得られる複数の測位衛星のそれぞれと受信機との距離から該受信機の位置を測位する単独測位方法であって、
測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値から測位衛星の軌道情報および衛星時計誤差を観測し、電離層遅延量情報、および対流圏遅延量情報を取得する対流圏遅延情報取得し、過去の受信機位置推定結果および測位衛星の軌道情報を用いて受信機位置を線形近似化し、未知数である線形近似された受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とする回帰方程式を構成し、
該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算することを特徴としている。

これらの構成では、搬送波毎で測位衛星毎のキャリア位相、ＰＮコード毎で測位衛星毎のコード擬似距離、測位衛星毎の衛星軌道情報、時計誤差、電離層遅延情報、対流圏遅延情報、を観測値とし、受信機の３次元位置、受信機の時計誤差、搬送波毎で測位衛星毎の整数値バイアスを未知数する。そして、各観測値を用いて目的変数を形成し、各未知数を用いて説明変数を形成することで回帰方程式が構成される。この際、受信機位置は、過去の受信機位置推定結果および測位衛星の軌道情報を用いて線形近似される。そして、この回帰方程式に最小二乗法等のパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を含む未知数が推定演算される。

また、この発明の単独測位装置および単独測位方法は、整数値バイアス推定法を用いて整数値バイアスをフィックスして受信機位置を推定演算することを特徴としている。

この構成では、搬送波の整数値バイアスに対してＬＡＭＢＤＡ法等の整数値バイアス推定法を適用することで、整数値バイアスが整数値としてフィックスされる。このフィックスされた整数値バイアスを回帰方程式に既知数として代入することで、回帰方程式内の未知数の数が減少し、推定精度が向上する。

また、この発明の単独測位装置および単独測位方法は、複数エポックで取得したキャリア位相、コード擬似距離、受信機時計誤差、衛星軌道情報、衛星時計誤差、電離層遅延量情報、および、対流圏遅延量情報を一時記憶手段で記憶し、測位演算手段で取得した複数エポックの各情報を用いて回帰方程式を構成することを特徴としている。

この構成では、複数エポックに亘る各情報を得ることで、観測値の数と未知数の数とが増加するが、観測値の数の増加の方が多いので、回帰方程式の未知数を推定がより容易になる。

また、この発明の単独測位装置および単独測位方法は、測位演算手段で、回帰方程式に状態推定アルゴリズムを適用することを特徴としている。

この構成では、前述の回帰方程式にカルマンフィルタや非線形フィルタ等の状態推定アルゴリズムを適用することで、移動する受信機位置の推定が可能となる。

また、この発明の単独測位装置および単独測位方法は、測位演算手段で、衛星軌道情報を目的変数に含み、衛星軌道誤差を説明変数に含む回帰方程式を用いて、測位衛星位置を推定演算することを特徴としている。

この構成では、測位衛星位置に対応する衛星軌道誤差を未知数として説明変数に加えることで、回帰方程式から測位衛星位置が推定演算される。

この発明によれば、電離層遅延量情報、対流圏遅延量情報、衛星軌道情報、衛星時計誤差を観測値として用い、且つキャリア位相およびコード擬似距離も観測値として用いて、受信機位置、整数値バイアスを未知数とする回帰方程式に適用することで、受信機位置および整数値バイアスの推定精度が向上する。これにより、１つの回帰方程式のみを用いるという簡素な演算で、受信機の単独測位を高精度に行うことができる。

また、この発明によれば、ＬＡＭＢＤＡ法等の整数値バイアス推定法を用いることで、より高精度に単独測位を行うことができる。

また、この発明によれば、複数エポックに亘る観測値を用いることで必要な観測衛星数を少なくして、高精度な単独測位を行うことができる。

また、この発明によれば、前述の線形回帰方程式にカルマンフィルタ等の状態推定アルゴリズムを適用することで、移動する受信機の位置を高精度に測位することができる。

また、この発明によれば、衛星軌道情報を観測値とし、衛星軌道誤差を未知数として推定演算することで、受信機と測位衛星との間の距離をより高精度に測位することができる。これにより、受信機位置をより高精度に測位することができる。

本発明の実施形態の単独測位装置の概略構成を示すブロック図である。 本発明の実施形態の単独測位装置とＧＰＳ受信機とからなる単独測位システムの処理フローを示すフローチャートである。 本発明の実施形態の単独測位方法を用いた場合の受信機位置（経度、緯度）を示した図である。 本発明の実施形態の単独測位方法を用いた場合の受信機位置（楕円体高）を示した図である。

符号の説明

１０−単独測位装置
１１−航法メッセージ解析部
１２−衛星情報処理部
１３−測位演算部
２０−ＧＰＳアンテナ
３０−ＧＰＳ受信機

本発明の実施形態に係る単独測位装置について図を参照して説明する。なお、以下の説明では、ＧＰＳについて説明するが、他の全てのＧＮＳＳ（全地球的航法衛星システム）に適用することができる。
図１は本実施形態の単独測位装置の概略構成を示すブロック図である。
また、図２は本実施形態の単独測位装置とＧＰＳ受信機とからなる単独測位システムの処理フローを示すフローチャートである。

図１に示すように、単独測位装置１０は、ＧＰＳ受信機３０に接続し、航法メッセージ解析部１１、衛星情報処理部１２、および測位演算部１３を備える。

ＧＰＳ受信機３０はアンテナ２０に接続し、アンテナ２０で受信したＧＰＳ衛星（測位衛星）からのＧＰＳ信号（測位用信号）より、既知の方法でＬ１波およびＬ２波のキャリア位相、Ｃ／ＡコードおよびＰ（Ｙ）コードのコード擬似距離（擬似距離）を取得するともに、Ｌ１波に重畳された航法メッセージを取得する（ｓ１→ｓ２）。そして、ＧＰＳ受信機３０は航法メッセージを航法メッセージ解析部１１に出力する。また、ＧＰＳ受信機３０はキャリア位相およびコード擬似距離を測位演算部１３に出力する。

航法メッセージ解析部１１は、入力された航法メッセージを解析して、電離層遅延量情報、各ＧＰＳ衛星の時計誤差、軌道情報を取得する。また、方法メッセージ解析部１１は、数式モデルを用いて対流圏遅延量情報を取得する。そして、航法メッセージ解析部１１は取得した各情報を衛星情報処理部１２に出力する。

衛星情報処理部１２は、ＧＰＳ衛星のエフェメリス情報を用いて測位に用いるＧＰＳ衛星を選定して、選定したＧＰＳ衛星に関する衛星軌道情報、衛星時計誤差情報、電離層遅延量情報、および対流圏量遅延情報を測位演算部１３に出力する。

測位演算部１３は、入力された各観測値、すなわち、選定された各ＧＰＳ衛星に対するＬ１波のキャリア位相、Ｌ２波のキャリア位相、Ｃ／Ａコードによるコード擬似距離、Ｐ（Ｙ）コードのコード擬似距離、各ＧＰＳ衛星の衛星軌道情報、各ＧＰＳ衛星の衛星時計誤差情報、電離層遅延量情報、対流圏遅延量情報を用いて、後述する線形回帰方程式を形成する。そして、測位演算部１３は、この線形回帰方程式にパラメータ推定アルゴリズムである最小二乗法を適用させて、Ｌ１、Ｌ２波に対する整数値バイアスＮ_L1，Ｎ_L2と受信機の位置とを推定演算する（ｓ４）。この演算は推定値の変化が予め設定された所定の閾値以下になるまで繰り返し行われ、推定値の変化が所定閾値に達した時点での受信機位置の推定演算結果が測位結果として出力される。

ここで、航法メッセージ解析部１１、衛星情報処理部１２、測位演算部１３は、それぞれ、以下に示すアルゴリズムを実現するマイクロプロセッサ等の数値演算処理器からなる。そして、これらの部分は複数の数値演算処理器により形成してもよく、１つの数値演算処理器により形成してもよい。

次に、前述の整数値バイアスＮ_L1，Ｎ_L2および受信機位置ｕの推定演算アルゴリズムについて詳述する。

一般に、受信機ｕ、ＧＰＳ衛星ｐに対するキャリア位相φ^p _L,uの観測方程式は式（１）で表すことができ、コード擬似距離（擬似距離）ρ^p _c,uの観測方程式は式（２）で表される。ここで、マルチパス誤差は微少として無視する。

ここで、λ_LはＬ波の波長を示す。また、ｒ^p _u（ｔ，ｔ−τ^p _u）は時刻ｔでの受信機ｕと時刻（ｔ−τ^p _u）でのＧＰＳ衛星との距離を示す。δＩ^p _u（ｔ）はＬ１波の電離層遅延量を示し、δＴ^p _u（ｔ）はＬ１波、Ｌ２波の対流圏量遅延を示す。δｔ_u（ｔ）は真の時刻ｔでの受信機ｕの時計誤差を示し、δｔ^p（ｔ−τ^p _u）は時刻（ｔ−τ^p _u）でのＧＰＳ衛星ｐの時計誤差を示す。Ｎ^p _L,uはＬ波における受信機ｕとＧＰＳ衛星ｐとの間の整数値バイアスを示し、ε^p _L,u（ｔ），ｅ^p _c,u（ｔ）はそれぞれ観測雑音を示す。

したがって、Ｌ１波におけるキャリア位相φ^p _L1,uの観測方程式は式（３）となり、Ｌ２波におけるキャリア位相φ^p _L2,uの観測方程式は式（４）となる。

ここで、ｆ_L1はＬ１波の周波数、ｆ_L2はＬ２の周波数を示す。

さらに、Φ^p _L1,u，Φ^p _L2,uをそれぞれ次式（３’），（４’）により定義する。

また、Ｃ／Ａコードによるコード擬似距離ρ^p _CA,uの観測方程式は式（５）となり、Ｐコードによるコード擬似距離ρ^p _P,uの観測方程式は式（６）となる。

ところで、受信機とＧＰＳ衛星との距離ｒ^p _u（ｔ，ｔ−τ^p _u）は、式（７）で表すことができる。

次に、未知数である受信機位置ｕ（ｔ）≡［ｘ_u（ｔ），ｙ_u（ｔ），ｚ_u（ｔ）］^Tを先験的な推定受信機位置ｕ^(j)（ｔ）≡［ｘ_u ^(j)（ｔ），ｙ_u ^(j)（ｔ），ｚ_u ^(j)（ｔ）］^Tのまわりで１次のティラー級数展開を行い、ｒ^p _u（ｔ）を線形近似すると次式を得られる。

この式（９）により、式（３’）、（４’）、（５）、（６）は、次式（１０）〜（１３）で表される。

ここで、Φe^p,(j) _L1,u，Φe^p,(j) _L2,u，ρe^p,(j) _L1,u，ρe^p,(j) _L2,u、を式（１４）〜（１７）で定義する。

これにより、式（１０）〜（１３）は次式（１８）〜（２１）となる。

これは、すなわち、キャリア位相、コード擬似距離を目的変数とし、受信機の位置、電離層遅延量、対流圏遅延量、整数値バイアス、誤差要因を説明変数とする近似的な線形回帰方程式に相当する。

ここで、Ｇ^(j)を式（２２）で定義する。

なお、Ｇ^(j) _uは式（２３）で表される。

また、Φe^p,(j) _L1,u，Φe^p,(j) _L2,u，ρe^p,(j) _L1,u，ρe^p,(j) _L2,u、を行列式で表す。

この式を用いることで、式（１８）〜式（２１）はベクトル行列からなる次式で表される。

もしくは、

ところで、前述のように、航法メッセージにはＧＰＳ衛星毎の時計誤差情報が含まれており、航法メッセージ解析部１１でＧＰＳ衛星毎の時計誤差情報が取得されている。この衛星時計誤差の観測値δｔｅ^pは次式で表され、前述の回帰方程式に加えることができる。

また、航法メッセージから得られるＧＰＳ衛星の情報からＧＰＳ衛星毎の電離層遅延量および対流圏遅延量をそれぞれ電離層遅延モデルおよび対流圏遅延モデルを用いて取得することができる。これにより、電離層遅延量の観測値δＩｅ_uとし、対流圏遅延量の観測値δＴｅ_uはそれぞれ次式で表され、前述の回帰方程式に加えることができる。

この結果、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量、対流圏遅延量を目的変数とし、整数値バイアス、受信機位置を説明変数として含む線形回帰方程式が実現される。

すなわち、式（２６）に示した回帰方程式は、次式に置き換えられる。

ここで、υ_uEの共分散行列をＲは式（３５）で表される。

従って、θ^(j+1)の推定値θｅ^(j+1)は、式（３６）となる。

また、推定値θｅ^(j+1)の分散値は、式（３７）となる。

この式（３６）に基づいて、最小二乗法を用いて繰り返し推定演算を行うことで、θを構成する未知数を収束させて各値を取得する。この際、説明変数ベクトルθを構成する未知数の収束条件は、未知数の推定値の繰り返し演算での変化量が、予め設定した閾値以下に達した時点である。そして、この条件下で説明変数ベクトルθの構成要素の値を確定させる。このような演算を行うことで、整数値バイアスＮ_L1，Ｎ_L2や受信機位置ｕを推定することができる。

このように、本実施形態の構成を用いることで、１重位相差等の処理を行うことなく整数値バイアスＮ_L1，Ｎ_L2や受信機位置ｕの推定演算、すなわち単独測位を行うことができる。これにより、従来例より簡素なアルゴリズムで単独測位を行うことができ、且つ単独測位の演算処理速度を向上することが可能になる。

ところで、このような推定演算を用いる場合、未知パラメータはθ_uの各要素であり、式（３１）に示す回帰方程式内の未知パラメータ数は、受信機位置のｕが３個、受信機時計誤差δｔ_uが１個、ＧＰＳ衛星の時計誤差δｔ^pがｎ_a個（ＧＰＳ衛星数）、電離層遅延δＩがｎ_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≧４を意味し、実質的に少なくとも４つの測位衛星からの測位用信号を受信すれば受信機の単独測位が可能となる。すなわち、少なくとも４つの測位衛星からの測位用信号を受信できれば、前述の単独測位は可能となる。

次に、本実施形態の単独測位方法を用いた場合のシミュレーション結果について説明する。なお、このシミュレーション結果は、Ｐコードを用いずに推定演算を行ったものである。

図３は本実施形態の単独測位方法を用いた場合の受信機位置（経度、緯度）を示した図である。図４は本実施形態の単独測位方法を用いた場合の受信機位置（楕円体高）を示した図である。

図３、図４において、円形プロットは本実施形態の単独測位結果１０１を示す。四角形プロットは従来のコード擬似距離のみを用いた単独測位結果１０２を示す。三角形プロットは従来のコード擬似距離、電離層遅延量、対流圏遅延量、および衛星時計誤差を用いた単独測位結果１０３を示す。また、図３におけるアスタリスクマーキング点および図４における実線は相対測位結果を示す。

また、表１は本実施形態の単独測位方法による受信機位置（経度、緯度、楕円体高）の平均および標準偏差を示す。表２は、従来のコード擬似距離のみを用いた単独測位方法による受信機位置（経度、緯度、楕円体高）の平均および標準偏差を示す。表３は、従来のコード擬似距離、電離層遅延量、対流圏遅延量、衛星時計誤差を用いた単独測位方法による受信機位置（経度、緯度、楕円体高）の平均および標準偏差を示す。

なお、このシミュレーションでは、受信機位置の初期値として、ＲＩＮＥＸデータの”ＡＰＰＲＯＸ ＰＯＳＩＴＩＯＮ ＸＹＺ”に記載された座標を用いた。また、電離層遅延量は、航法メッセージに含まれる電離層遅延に関するデータから、いわゆる放送モデル（Ｋｌｏｂｕｃｈａｒモデル）を適用して算出した値を用い、対流圏遅延量は、航法データから計算される測位衛星の仰角情報から次式を適用して算出した値を用いている。

ここで、ξは仰角である。
さらに、観測誤差の分散は、それぞれ、コード擬似距離が１．５［ｍ］（１σ）で、キャリア位相が（波長λ／１０）＋（１．５／１０）［ｍ］（１σ）で、衛星時計誤差が３．６［ｍ］（１σ）で、電離層遅延が７．０［ｍ］（１σ）で、対流圏遅延が０．７［ｍ］（１σ）でそれぞれ設定されている。また、このシミュレーションは、最小二乗法による繰り返し演算での、受信機位置の推定値の変化のノルム値が１×１０^-3［ｍ］以下になるまで計算を繰り返した。

これらの結果に示すように、本実施形態の単独測位方法を用いることにより、高精度に受信機位置を測位することができる。

以上のように、本実施形態の単独測位装置および単独測位方法を用いることにより、従来よりも簡素な演算処理で高精度に受信機の単独測位を行うことができる。

なお、前述の最小二乗法による推定演算時に、式（３７）より得られる分散値を用いて、ＬＡＭＢＤＡ法を適用することで、整数値バイアスＮ_L1，Ｎ_L2をフィックスして、さらに高精度の受信機位置を測位することができる。

また、前述の説明では、１エポックの観測値を用いて単独測位を行う方法を示したが複数エポックの観測値を用いても単独測位を行うことができる。この際、単独測位装置には、複数エポックのキャリア位相、コード擬似距離、測位衛星の衛星軌道情報、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を記憶するメモリ（本発明の「一時記憶手段」に相当する。）を備える。

このような場合、例えば、Ｃ／Ａコードのみを使用した場合、Ｋエポックで取得できる観測値は、次式で表される。

６ｎ_s×Ｋ＝６Ｋｎ_s
一方、観測エポックが複数となるキネマティック測位では、整数値バイアスＮ_L1，Ｎ_L2を除く各未知数（未知パラメータ）がエポック数に応じて増加するので、未知数の数は、次式で表される。

４＋５ｎ_s＋（Ｋ−１）（４＋３ｎ_s）＝４Ｋ＋（３Ｋ＋２）ｎ_s
したがって、回帰方程式を最小二乗法により解くことができる必要条件は、次式となる。

６Ｋｎ_s≧４Ｋ＋（３Ｋ＋２）ｎ_s
（３Ｋ−２）ｎ_s≧４Ｋ
すなわち、
ｎ_s≧４Ｋ／（３Ｋ−２）
そして、エポック数Ｋを２→∞の場合、ｎ_sの取り得る整数値を計算すると、ＧＰＳ衛星数ｎ_sは次式の値となる。

ｎ_s≧２
これにより、観測エポック数を増加させることにより、観測するＧＰＳ衛星数を最低２個として測位を行なうことができる。

ここで、このようなキネマティック測位において、カルマンフィルタを適用することで、より高精度な測位が可能となる。以下に、カルマンフィルタを適用した場合の測位アルゴリズムについて説明する。

まず、式（３１）より、エポック数に対応する観測時間ｔを要素として表現し直すと、

これを、ベクトル行列表示すると、式（３９）になる。

ここでｃδｔ^p，δＩ_u，δＴ_uをそれぞれ式（４０ａ），（４０ｂ），（４０ｃ）とする。

これらを式（３９）の回帰方程式に代入し、さらに、次式に示す定義を行う。

この結果、新たな回帰方程式は、式（４５）で表される。

ここで、未知の受信機位置ｕ（ｔ）の速度に対するマルコフ過程モデルを仮定し、かつ、受信機の時計誤差ｃδｔ_u（ｔ）に対してマルコフ過程モデルを仮定できるものとし、速度ベクトルをｖ（ｔ）とおき、受信機の新たな状態ベクトルη_u（ｔ）を次式で定義する。

そして、この状態ベクトルη_u（ｔ）について次の状態方程式を構成する。

この際、式（４５）の観測方程式は、次式で表すことができる。

この観測方程式をベクトル表示すると、式（４９）となる。

これは、すなわち、式（４７）、式（４８）（式（４９））からなるカルマンフィルタを意味する。この際、観測誤差υ_u,R（ｔ）の共分散行列Ｒは、式（３５）と同様に求めることができる。そして、これら状態方程式（４７）、観測方程式（４８）、を用いることにより移動する受信機位置を推定演算することができる。

このような構成とすることで、受信機位置が移動してもカルマンフィルタにより受信機の移動を推定しながら、受信機位置を高精度に測位することができる。

ところで、前述の説明では、ＧＰＳ衛星の位置を推定演算しなかったが、次に示す方法（アルゴリズム）を用いることにより、ＧＰＳ衛星の位置を推定演算することができる。

式（７）に示すように、受信機とＧＰＳ衛星との距離を定義し、受信機とともにＧＰＳ衛星の位置を推定演算する場合、受信機位置に関する線形近似式は、前述のように表される。

測位衛星位置に関する線形近似式は、次式で表される。

そして、ｓｐ≡［ｘ^p，ｙ^p，ｚ^p］^Tと定義すると次式の関係が成り立つ。

したがって、ｕ，ｓｐの推定値は、それぞれｕ^(j)，ｓｅ^pと１次ティラー級数展開とを用いて次式で近似される。

したがって、式（１０）から式（１３）は、それぞれ次式で表される。

ここで、Φe^p,(j) _L1,u，Φe^p,(j) _L2,u，ρe^p,(j) _L1,u，ρe^p,(j) _L2,u、をそれぞれ式（５５）〜（５８）で定義する。

これにより、式（５１）〜（５４）は次式（５９）〜（６２）となる。

ここで、Φe^p,(j) _L1,u，Φe^p,(j) _L2,u，ρe^p,(j) _L1,u，ρe^p,(j) _L2,u、をそれぞれ次式で定義する。

これにより、式（５９）〜式（６２）はベクトル行列からなる次式で表される。

また、精密軌道衛星に関する観測値ベクトルｓｋを次式で表す。

これを式（６３）に適用すると式（６４）となる。

このような回帰方程式を用いることで、ＧＰＳ衛星の位置についても推定演算することができる。そして、ＧＰＳ衛星位置の推定位置精度が向上することで、受信機位置の推定位置をより高精度に推定することができる。

さらに、以下のアルゴズムを用いることで、このような場合でもカルマンフィルタを適用することができる。

すなわち、ｓ，ｃδｔ^p，δＩ_u，δＴ_uをそれぞれ式（６５ａ）〜（６５ｄ）で定義する。

これらの関係を用い、式（６４）の回帰方程式に代入し、さらに
次に示す定義を行う。

この結果、新たな回帰方程式は、式（７０）で表される。

そして、この観測方程式に対して前述のように状態方程式を定義することで、カルマンフィルタを構成することができる。

なお、前述の実施形態では、受信機で電離層遅延量や対流圏遅延量等を推定していたが、位置が既知であり固定されている基地局が存在すれば、基地局において、前述の方法を用いて電離層遅延量や対流圏遅延量を推定しても良い。そして、このように基地局で推定した電離層遅延量や対流圏遅延量、さらには、衛星軌道誤差、衛星時計誤差、を受信機に与えることにより、受信機でさらに高精度な推定を行うことができる。

また、前述の実施形態では、線形回帰方程式に最小二乗法を適用した例を示したが、他のパラメータ推定アルゴリズムを用いてもよい。

また、前述の実施形態では、受信機位置の線形近似に１次のティラー級数展開を用いたが、他の線形化演算を用いてもよい。

また、前述の実施形態では、ＬＡＭＢＤＡ法を利用する例を示したが、他の整数値バイアス推定方法を用いてもよい。

また、前述の実施形態では、カルマンフィルタを適用する例を示したが、他の状態推定アルゴリズムを用いてもよい。

本発明は、測位用信号を用いて受信機の単独測位を行う単独測位装置および単独測位方法に関するものである。

従来、測位衛星から送信される測位用信号を用いて単独測位を行う装置が各種開示されている。これらの装置による基本測位演算は、コード擬似距離を用いて、受信機の３次元の位置（誤差）と受信機の時計誤差とを未知数とした非線形連立方程式にニュートン法や拡張カルマンフィルタを適用することで求めるものであった。この際、電離層遅延による誤差や対流圏遅延による誤差の影響を除去する方法が各種考案されている。その１つとして、これらの遅延による影響を初期状態から「０」に設定する方法があり、またその１つとして、１重位相差等の演算処理を行ってこれらの遅延量を推定する方法があった。
日本測地学会編著，「新訂版ＧＰＳ−人工衛星による精密測位システム−」，社団法人日本測量協会，１９８９年１１月１５日，ｐ．１２１−１４０

ところが、従来のコード擬似距離のみを用いた単独測位装置および単独測位方法では、現実に存在する電離層遅延や対流圏遅延を無視したり、あるいは十分な精度で求めずして推定演算を行っていた。このため、測位の推定値のばらつきが大きく、高精度な測位を行うことができなかった。また、コード擬似距離を用いるとともに、電離層遅延情報や対流圏遅延情報を用いて１重位相差処理等を行う場合、前述のコード擬似距離のみを用いる場合よりも推定精度は向上するものの、推定演算内の１重位相差処理を行う等、推定演算処理が複雑になったり、１重位相差処理を行った後のノイズの設定により大きく推定演算結果が異なったりする。このため、推定演算の処理内容が複雑なわりには、推定演算結果の精度の向上性が悪く、さらに推定演算処理速度が遅くなる可能性があった。

したがって、この発明の目的は、推定演算処理を複雑にすることなく確実に高精度測位を実現することができる単独測位装置および単独測位方法を提供することにある。

この発明は、複数の測位衛星から送信される測位用信号を用いて得られる複数の測位衛星のそれぞれと受信機との距離から該受信機の位置を測位する単独測位装置において、
測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値から測位衛星の軌道情報および衛星時計誤差を観測する衛星情報検出手段と、電離層遅延量情報を取得する電離層遅延情報取得手段と、対流圏遅延量情報を取得する対流圏遅延情報取得手段と、受信機位置を過去の受信機位置推定結果および測位衛星の軌道情報を用いて線形近似化し、未知数である線形近似された受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とする回帰方程式を構成し、該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算する測位演算手段と、を備えたことを特徴としている。

また、この発明は、複数の測位衛星から送信される測位用信号を用いて得られる複数の測位衛星のそれぞれと受信機との距離から該受信機の位置を測位する単独測位方法であって、
測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値から測位衛星の軌道情報および衛星時計誤差を観測し、電離層遅延量情報、および対流圏遅延量情報を取得する対流圏遅延情報取得し、過去の受信機位置推定結果および測位衛星の軌道情報を用いて受信機位置を線形近似化し、未知数である線形近似された受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とする回帰方程式を構成し、
該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算することを特徴としている。

これらの構成では、搬送波毎で測位衛星毎のキャリア位相、ＰＮコード毎で測位衛星毎のコード擬似距離、測位衛星毎の衛星軌道情報、時計誤差、電離層遅延情報、対流圏遅延情報、を観測値とし、受信機の３次元位置、受信機の時計誤差、搬送波毎で測位衛星毎の整数値バイアスを未知数する。そして、各観測値を用いて目的変数を形成し、各未知数を用いて説明変数を形成することで回帰方程式が構成される。この際、受信機位置は、過去の受信機位置推定結果および測位衛星の軌道情報を用いて線形近似される。そして、この回帰方程式に最小二乗法等のパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を含む未知数が推定演算される。

また、この発明の単独測位装置および単独測位方法は、整数値バイアス推定法を用いて整数値バイアスをフィックスして受信機位置を推定演算することを特徴としている。

この構成では、搬送波の整数値バイアスに対してＬＡＭＢＤＡ法等の整数値バイアス推定法を適用することで、整数値バイアスが整数値としてフィックスされる。このフィックスされた整数値バイアスを回帰方程式に既知数として代入することで、回帰方程式内の未知数の数が減少し、推定精度が向上する。

また、この発明の単独測位装置および単独測位方法は、複数エポックで取得したキャリア位相、コード擬似距離、受信機時計誤差、衛星軌道情報、衛星時計誤差、電離層遅延量情報、および、対流圏遅延量情報を一時記憶手段で記憶し、測位演算手段で取得した複数エポックの各情報を用いて回帰方程式を構成することを特徴としている。

この構成では、複数エポックに亘る各情報を得ることで、観測値の数と未知数の数とが増加するが、観測値の数の増加の方が多いので、回帰方程式の未知数を推定がより容易になる。

また、この発明の単独測位装置および単独測位方法は、測位演算手段で、回帰方程式に状態推定アルゴリズムを適用することを特徴としている。

この構成では、前述の回帰方程式にカルマンフィルタや非線形フィルタ等の状態推定アルゴリズムを適用することで、移動する受信機位置の推定が可能となる。

また、この発明の単独測位装置および単独測位方法は、測位演算手段で、衛星軌道情報を目的変数に含み、衛星軌道誤差を説明変数に含む回帰方程式を用いて、測位衛星位置を推定演算することを特徴としている。

この構成では、測位衛星位置に対応する衛星軌道誤差を未知数として説明変数に加えることで、回帰方程式から測位衛星位置が推定演算される。

この発明によれば、電離層遅延量情報、対流圏遅延量情報、衛星軌道情報、衛星時計誤差を観測値として用い、且つキャリア位相およびコード擬似距離も観測値として用いて、受信機位置、整数値バイアスを未知数とする回帰方程式に適用することで、受信機位置および整数値バイアスの推定精度が向上する。これにより、１つの回帰方程式のみを用いるという簡素な演算で、受信機の単独測位を高精度に行うことができる。

また、この発明によれば、ＬＡＭＢＤＡ法等の整数値バイアス推定法を用いることで、より高精度に単独測位を行うことができる。

また、この発明によれば、複数エポックに亘る観測値を用いることで必要な観測衛星数を少なくして、高精度な単独測位を行うことができる。

また、この発明によれば、前述の線形回帰方程式にカルマンフィルタ等の状態推定アルゴリズムを適用することで、移動する受信機の位置を高精度に測位することができる。

また、この発明によれば、衛星軌道情報を観測値とし、衛星軌道誤差を未知数として推定演算することで、受信機と測位衛星との間の距離をより高精度に測位することができる。これにより、受信機位置をより高精度に測位することができる。

本発明の実施形態に係る単独測位装置について図を参照して説明する。なお、以下の説明では、ＧＰＳについて説明するが、他の全てのＧＮＳＳ（全地球的航法衛星システム）に適用することができる。
図１は本実施形態の単独測位装置の概略構成を示すブロック図である。
また、図２は本実施形態の単独測位装置とＧＰＳ受信機とからなる単独測位システムの処理フローを示すフローチャートである。

図１に示すように、単独測位装置１０は、ＧＰＳ受信機３０に接続し、航法メッセージ解析部１１、衛星情報処理部１２、および測位演算部１３を備える。

ＧＰＳ受信機３０はアンテナ２０に接続し、アンテナ２０で受信したＧＰＳ衛星（測位衛星）からのＧＰＳ信号（測位用信号）より、既知の方法でＬ１波およびＬ２波のキャリア位相、Ｃ／ＡコードおよびＰ（Ｙ）コードのコード擬似距離（擬似距離）を取得するともに、Ｌ１波に重畳された航法メッセージを取得する（ｓ１→ｓ２）。そして、ＧＰＳ受信機３０は航法メッセージを航法メッセージ解析部１１に出力する。また、ＧＰＳ受信機３０はキャリア位相およびコード擬似距離を測位演算部１３に出力する。

航法メッセージ解析部１１は、入力された航法メッセージを解析して、電離層遅延量情報、各ＧＰＳ衛星の時計誤差、軌道情報を取得する。また、方法メッセージ解析部１１は、数式モデルを用いて対流圏遅延量情報を取得する。そして、航法メッセージ解析部１１は取得した各情報を衛星情報処理部１２に出力する。

衛星情報処理部１２は、ＧＰＳ衛星のエフェメリス情報を用いて測位に用いるＧＰＳ衛星を選定して、選定したＧＰＳ衛星に関する衛星軌道情報、衛星時計誤差情報、電離層遅延量情報、および対流圏量遅延情報を測位演算部１３に出力する。

測位演算部１３は、入力された各観測値、すなわち、選定された各ＧＰＳ衛星に対するＬ１波のキャリア位相、Ｌ２波のキャリア位相、Ｃ／Ａコードによるコード擬似距離、Ｐ（Ｙ）コードのコード擬似距離、各ＧＰＳ衛星の衛星軌道情報、各ＧＰＳ衛星の衛星時計誤差情報、電離層遅延量情報、対流圏遅延量情報を用いて、後述する線形回帰方程式を形成する。そして、測位演算部１３は、この線形回帰方程式にパラメータ推定アルゴリズムである最小二乗法を適用させて、Ｌ１、Ｌ２波に対する整数値バイアスＮL1，ＮL2と受信機の位置とを推定演算する（ｓ４）。この演算は推定値の変化が予め設定された所定の閾値以下になるまで繰り返し行われ、推定値の変化が所定閾値に達した時点での受信機位置の推定演算結果が測位結果として出力される。

ここで、航法メッセージ解析部１１、衛星情報処理部１２、測位演算部１３は、それぞれ、以下に示すアルゴリズムを実現するマイクロプロセッサ等の数値演算処理器からなる。そして、これらの部分は複数の数値演算処理器により形成してもよく、１つの数値演算処理器により形成してもよい。

次に、前述の整数値バイアスＮ_L1，Ｎ_L2および受信機位置ｕの推定演算アルゴリズムについて詳述する。

一般に、受信機ｕ、ＧＰＳ衛星ｐに対するキャリア位相φ^p _L,uの観測方程式は式（１）で表すことができ、コード擬似距離（擬似距離）ρ^p _c,uの観測方程式は式（２）で表される。ここで、マルチパス誤差は微少として無視する。

ここで、λLはＬ波の波長を示し、ｒ^p _u（ｔ，ｔ−τ^p _u）は時刻ｔでの受信機ｕと時刻（ｔ−τ^p _u）でのＧＰＳ衛星との距離を示す。δＩ^p _u（ｔ）はＬ１波の電離層遅延量を示し、δＴ^p _u（ｔ）はＬ１波、Ｌ２波の対流圏量遅延を示す。δｔ_u（ｔ）は真の時刻ｔでの受信機ｕの時計誤差を示し、δｔ^p（ｔ−τ^p _u）は時刻（ｔ−τ^p _u）でのＧＰＳ衛星ｐの時計誤差を示す。Ｎ^p _L,uはＬ波における受信機ｕとＧＰＳ衛星ｐとの間の整数値バイアスを示し、ε^p _L,u（ｔ），ｅ^p _c,u（ｔ）はそれぞれ観測雑音を示す。

したがって、Ｌ１波におけるキャリア位相φ^p _L1,uの観測方程式は式（３）となり、Ｌ２波におけるキャリア位相φ^p _L2,uの観測方程式は式（４）となる。

ここで、ｆ_L1はＬ１波の周波数、ｆ_L2はＬ２の周波数を示す。

さらに、Φ^p _L1,u，Φ^p _L2,uをそれぞれ次式（３’），（４’）により定義する。

また、Ｃ／Ａコードによるコード擬似距離ρ^p _CA,uの観測方程式は式（５）となり、Ｐコードによるコード擬似距離ρ^p _P,uの観測方程式は式（６）となる。

ところで、受信機とＧＰＳ衛星との距離ｒ^p _u（ｔ，ｔ−τ^p _u）は、式（７）で表すことができる。

次に、未知数である受信機位置ｕ（ｔ）≡［ｘ_u（ｔ），ｙ_u（ｔ），ｚ_u（ｔ）］Tを先験的な推定受信機位置ｕ^(j)（ｔ）≡［ｘ_u ^(j)（ｔ），ｙ_u ^(j)（ｔ），ｚ_u ^(j)（ｔ）］Tのまわりで１次のティラー級数展開を行い、ｒ^p _u（ｔ）を線形近似すると次式を得られる。

この式（９）により、式（３’）、（４’）、（５）、（６）は、次式（１０）〜（１３）で表される。

ここで、Φe^p,(j) _L1,u，Φe^p,(j) _L2,u，ρe^p,(j) _L1,u，ρe^p,(j) _L2,u、を式（１４）〜（１７）で定義する。

これにより、式（１０）〜（１３）は次式（１８）〜（２１）となる。

これは、すなわち、キャリア位相、コード擬似距離を目的変数とし、受信機の位置、電離層遅延量、対流圏遅延量、整数値バイアス、誤差要因を説明変数とする近似的な線形回帰方程式に相当する。

ここで、Ｇ^(j)を式（２２）で定義する。

なお、Ｇ^(j) _uは式（２３）で表される。

また、Φe^p,(j) _L1,u，Φe^p,(j) _L2,u，ρe^p,(j) _L1,u，ρe^p,(j) _L2,u、を行列式で表す。

この式を用いることで、式（１８）〜式（２１）はベクトル行列からなる次式で表される。

もしくは、

ところで、前述のように、航法メッセージにはＧＰＳ衛星毎の時計誤差情報が含まれており、航法メッセージ解析部１１でＧＰＳ衛星毎の時計誤差情報が取得されている。この衛星時計誤差の観測値δｔｅ^pは次式で表され、前述の回帰方程式に加えることができる。

また、航法メッセージから得られるＧＰＳ衛星の情報からＧＰＳ衛星毎の電離層遅延量および対流圏遅延量をそれぞれ電離層遅延モデルおよび対流圏遅延モデルを用いて取得することができる。これにより、電離層遅延量の観測値δＩｅ_uとし、対流圏遅延量の観測値δＴｅ_uはそれぞれ次式で表され、前述の回帰方程式に加えることができる。

この結果、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量、対流圏遅延量を目的変数とし、整数値バイアス、受信機位置を説明変数として含む線形回帰方程式が実現される。

すなわち、式（２６）に示した回帰方程式は、次式に置き換えられる。

ここで、υ_uEの共分散行列をＲは式（３５）で表される。

従って、θ^(j+1)の推定値θｅ^(j+1)は、式（３６）となる。

また、推定値θｅ^(j+1)の分散値は、式（３７）となる。

この式（３６）に基づいて、最小二乗法を用いて繰り返し推定演算を行うことで、θを構成する未知数を収束させて各値を取得する。この際、説明変数ベクトルθを構成する未知数の収束条件は、未知数の推定値の繰り返し演算での変化量が、予め設定した閾値以下に達した時点である。そして、この条件下で説明変数ベクトルθの構成要素の値を確定させる。このような演算を行うことで、整数値バイアスＮ_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≧４を意味し、実質的に少なくとも４つの測位衛星からの測位用信号を受信すれば受信機の単独測位が可能となる。すなわち、少なくとも４つの測位衛星からの測位用信号を受信できれば、前述の単独測位は可能となる。

次に、本実施形態の単独測位方法を用いた場合のシミュレーション結果について説明する。なお、このシミュレーション結果は、Ｐコードを用いずに推定演算を行ったものである。

図３は本実施形態の単独測位方法を用いた場合の受信機位置（経度、緯度）を示した図である。図４は本実施形態の単独測位方法を用いた場合の受信機位置（楕円体高）を示した図である。

図３、図４において、円形プロットは本実施形態の単独測位結果１０１を示す。四角形プロットは従来のコード擬似距離のみを用いた単独測位結果１０２を示す。三角形プロットは従来のコード擬似距離、電離層遅延量、対流圏遅延量、および衛星時計誤差を用いた単独測位結果１０３を示す。また、図３におけるアスタリスクマーキング点および図４における実線は相対測位結果を示す。

また、表１は本実施形態の単独測位方法による受信機位置（経度、緯度、楕円体高）の平均および標準偏差を示す。表２は、従来のコード擬似距離のみを用いた単独測位方法による受信機位置（経度、緯度、楕円体高）の平均および標準偏差を示す。表３は、従来のコード擬似距離、電離層遅延量、対流圏遅延量、衛星時計誤差を用いた単独測位方法による受信機位置（経度、緯度、楕円体高）の平均および標準偏差を示す。

なお、このシミュレーションでは、受信機位置の初期値として、ＲＩＮＥＸデータの”ＡＰＰＲＯＸ ＰＯＳＩＴＩＯＮ ＸＹＺ”に記載された座標を用いた。また、電離層遅延量は、航法メッセージに含まれる電離層遅延に関するデータから、いわゆる放送モデル（Ｋｌｏｂｕｃｈａｒモデル）を適用して算出した値を用い、対流圏遅延量は、航法データから計算される測位衛星の仰角情報から次式を適用して算出した値を用いている。

ここで、ξは仰角である。
さらに、観測誤差の分散は、それぞれ、コード擬似距離が１．５［ｍ］（１σ）で、キャリア位相が（波長λ／１０）＋（１．５／１０）［ｍ］（１σ）で、衛星時計誤差が３．６［ｍ］（１σ）で、電離層遅延が７．０［ｍ］（１σ）で、対流圏遅延が０．７［ｍ］（１σ）でそれぞれ設定されている。また、このシミュレーションは、最小二乗法による繰り返し演算での、受信機位置の推定値の変化のノルム値が１×１０^-3［ｍ］以下になるまで計算を繰り返した。

これらの結果に示すように、本実施形態の単独測位方法を用いることにより、高精度に受信機位置を測位することができる。

以上のように、本実施形態の単独測位装置および単独測位方法を用いることにより、従来よりも簡素な演算処理で高精度に受信機の単独測位を行うことができる。

なお、前述の最小二乗法による推定演算時に、式（３７）より得られる分散値を用いて、ＬＡＭＢＤＡ法を適用することで、整数値バイアスＮ_L1，Ｎ_L2をフィックスして、さらに高精度の受信機位置を測位することができる。

また、前述の説明では、１エポックの観測値を用いて単独測位を行う方法を示したが複数エポックの観測値を用いても単独測位を行うことができる。この際、単独測位装置には、複数エポックのキャリア位相、コード擬似距離、測位衛星の衛星軌道情報、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を記憶するメモリ（本発明の「一時記憶手段」に相当する。）を備える。

このような場合、例えば、Ｃ／Ａコードのみを使用した場合、Ｋエポックで取得できる観測値は、次式で表される。

６ｎ_s×Ｋ＝６Ｋｎ_s
一方、観測エポックが複数となるキネマティック測位では、整数値バイアスＮ_L1，Ｎ_L2を除く各未知数（未知パラメータ）がエポック数に応じて増加するので、未知数の数は、次式で表される。

４＋５ｎ_s＋（Ｋ−１）（４＋３ｎ_s）＝４Ｋ＋（３Ｋ＋２）ｎ_s
したがって、回帰方程式を最小二乗法により解くことができる必要条件は、次式となる。

６Ｋｎ_s≧４Ｋ＋（３Ｋ＋２）ｎ_s
（３Ｋ−２）ｎ_s≧４Ｋ
すなわち、
ｎ_s≧４Ｋ／（３Ｋ−２）
そして、エポック数Ｋを２→∞の場合、ｎ_sの取り得る整数値を計算すると、ＧＰＳ衛星数ｎ_sは次式の値となる。

ｎ_s≧２
これにより、観測エポック数を増加させることにより、観測するＧＰＳ衛星数を最低２個として測位を行なうことができる。

ここで、このようなキネマティック測位において、カルマンフィルタを適用することで、より高精度な測位が可能となる。以下に、カルマンフィルタを適用した場合の測位アルゴリズムについて説明する。

まず、式（３１）より、エポック数に対応する観測時間ｔを要素として表現し直すと、

これを、ベクトル行列表示すると、式（３９）になる。

ここでｃδｔ^p，δＩ_u，δＴ_uをそれぞれ式（４０ａ），（４０ｂ），（４０ｃ）とする。

これらを式（３９）の回帰方程式に代入し、さらに、次式に示す定義を行う。

この結果、新たな回帰方程式は、式（４５）で表される。

ここで、未知の受信機位置ｕ（ｔ）の速度に対するマルコフ過程モデルを仮定し、かつ、受信機の時計誤差ｃδｔ_u（ｔ）に対してマルコフ過程モデルを仮定できるものとし、速度ベクトルをｖ（ｔ）とおき、受信機の新たな状態ベクトルη_u（ｔ）を次式で定義する。

そして、この状態ベクトルη_u（ｔ）について次の状態方程式を構成する。

この際、式（４５）の観測方程式は、次式で表すことができる。

この観測方程式をベクトル表示すると、式（４９）となる。

これは、すなわち、式（４７）、式（４８）（式（４９））からなるカルマンフィルタを意味する。この際、観測誤差υ_u,R（ｔ）の共分散行列Ｒは、式（３５）と同様に求めることができる。そして、これら状態方程式（４７）、観測方程式（４８）、を用いることにより移動する受信機位置を推定演算することができる。

このような構成とすることで、受信機位置が移動してもカルマンフィルタにより受信機の移動を推定しながら、受信機位置を高精度に測位することができる。

ところで、前述の説明では、ＧＰＳ衛星の位置を推定演算しなかったが、次に示す方法（アルゴリズム）を用いることにより、ＧＰＳ衛星の位置を推定演算することができる。

式（７）に示すように、受信機とＧＰＳ衛星との距離を定義し、受信機とともにＧＰＳ衛星の位置を推定演算する場合、受信機位置に関する線形近似式は、前述のように表される。

測位衛星位置に関する線形近似式は、次式で表される。

そして、ｓｐ≡［ｘ^p，ｙ^p，ｚ^p］Tと定義すると次式の関係が成り立つ。

したがって、ｕ，ｓｐの推定値は、それぞれｕ^(j)，ｓｅ^pと１次ティラー級数展開とを用いて次式で近似される。

したがって、式（１０）から式（１３）は、それぞれ次式で表される。

ここで、Φe^p,(j) _L1,u，Φe^p,(j) _L2,u，ρe^p,(j) _L1,u，ρe^p,(j) _L2,u、をそれぞれ式（５５）〜（５８）で定義する。

これにより、式（５１）〜（５４）は次式（５９）〜（６２）となる。

ここで、Φe^p,(j) _L1,u，Φe^p,(j) _L2,u，ρe^p,(j) _L1,u，ρe^p,(j) _L2,u、をそれぞれ次式で定義する。

これにより、式（５９）〜式（６２）はベクトル行列からなる次式で表される。

また、精密軌道衛星に関する観測値ベクトルｓｋを次式で表す。

これを式（６３）に適用すると式（６４）となる。

このような回帰方程式を用いることで、ＧＰＳ衛星の位置についても推定演算することができる。そして、ＧＰＳ衛星位置の推定位置精度が向上することで、受信機位置の推定位置をより高精度に推定することができる。

さらに、以下のアルゴズムを用いることで、このような場合でもカルマンフィルタを適用することができる。

すなわち、ｓ，ｃδｔ^p，δＩ_u，δＴ_uをそれぞれ式（６５ａ）〜（６５ｄ）で定義する。

これらの関係を用い、式（６４）の回帰方程式に代入し、さらに
次に示す定義を行う。

この結果、新たな回帰方程式は、式（７０）で表される。

そして、この観測方程式に対して前述のように状態方程式を定義することで、カルマンフィルタを構成することができる。

なお、前述の実施形態では、線形回帰方程式に最小二乗法を適用した例を示したが、他のパラメータ推定アルゴリズムを用いてもよい。

また、前述の実施形態では、受信機位置の線形近似に１次のティラー級数展開を用いたが、他の線形化演算を用いてもよい。

また、前述の実施形態では、ＬＡＭＢＤＡ法を利用する例を示したが、他の整数値バイアス推定方法を用いてもよい。

また、前述の実施形態では、カルマンフィルタを適用する例を示したが、他の状態推定アルゴリズムを用いてもよい。

本発明の実施形態の単独測位装置の概略構成を示すブロック図である。 本発明の実施形態の単独測位装置とＧＰＳ受信機とからなる単独測位システムの処理フローを示すフローチャートである。 本発明の実施形態の単独測位方法を用いた場合の受信機位置（経度、緯度）を示した図である。 本発明の実施形態の単独測位方法を用いた場合の受信機位置（楕円体高）を示した図である。

符号の説明

１０−単独測位装置
１１−航法メッセージ解析部
１２−衛星情報処理部
１３−測位演算部
２０−ＧＰＳアンテナ
３０−ＧＰＳ受信機

Claims (10)

1. 複数の測位衛星から送信される測位用信号を用いて得られる前記複数の測位衛星のそれぞれと受信機との距離から該受信機の位置を測位する単独測位装置において、
   前記測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値から前記測位衛星の軌道情報および衛星時計誤差を観測する衛星情報検出手段と、
   電離層遅延量情報を取得する電離層遅延情報取得手段と、
   対流圏遅延量情報を取得する対流圏遅延情報取得手段と、
   受信機位置を過去の受信機位置推定結果および前記測位衛星の軌道情報を用いて線形近似化し、
   未知数である線形近似された受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とする回帰方程式を構成し、該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算する測位演算手段と、
   を備えたことを特徴とする単独測位装置。
2. 前記測位演算手段は、整数値バイアス推定法を用いて前記整数値バイアスをフィックスして、前記受信機位置を推定演算する請求項１に記載の単独測位装置。
3. 複数エポックで取得した前記キャリア位相、前記コード擬似距離、前記受信機時計誤差、前記衛星軌道情報、前記衛星時計誤差、前記電離層遅延量情報、および、前記対流圏遅延量情報を記憶する一時記憶手段を備え、
   前記測位演算手段は、取得した前記複数エポックの各情報を用いて、前記回帰方程式を構成する請求項１または請求項２に記載の単独測位装置。
4. 前記測位演算手段は、前記回帰方程式に状態推定アルゴリズムを適用する請求項１〜３のいずれかに記載の単独測位装置。
5. 前記測位演算手段は、前記衛星軌道情報を目的変数に含み、前記衛星軌道情報に対応する衛星軌道誤差を説明変数に含む前記回帰方程式を用いて、測位衛星位置を推定演算する請求項１〜４のいずれかに記載の単独測位装置。
6. 複数の測位衛星から送信される測位用信号を用いて得られる前記複数の測位衛星のそれぞれと受信機との距離から該受信機の位置を測位する単独測位方法であって、
   前記測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値から前記測位衛星の軌道情報および衛星時計誤差を観測し、
   電離層遅延量情報、および対流圏遅延量情報を取得する対流圏遅延情報取得し、
   過去の受信機位置推定結果および前記測位衛星の軌道情報を用いて受信機位置を線形近似化し、
   未知数である線形近似された受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とする回帰方程式を構成し、
   該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算する、
   ことを特徴とする単独測位方法。
7. 整数値バイアス推定法を用いて前記整数値バイアスをフィックスして、前記受信機位置を推定演算する請求項６に記載の単独測位方法。
8. 複数エポックで取得した前記キャリア位相、前記コード擬似距離、前記受信機時計誤差、前記衛星軌道情報、前記衛星時計誤差、前記電離層遅延量情報、および、前記対流圏遅延量情報を一時記憶し、
   取得した前記複数エポックの各情報を用いて、前記回帰方程式を構成する請求項６または請求項７に記載の単独測位方法。
9. 前記回帰方程式に状態推定アルゴリズムを適用する請求項６〜８のいずれかに記載の単独測位方法。
10. 前記衛星軌道情報を目的変数に含み、前記衛星軌道情報に対応する衛星軌道誤差を説明変数に含む前記回帰方程式を用いて、測位衛星位置を推定演算する請求項６〜９のいずれかに記載の単独測位方法。

Images