WO2006022318A1

WO2006022318A1 - 単独測位装置および単独測位方法

Info

Publication number: WO2006022318A1
Application number: PCT/JP2005/015404
Authority: WO
Inventors: Sueo Sugimoto; Yukihiro Kubo; Kazunori Kagawa; Yasuhiro Tajima; Mitsuru Nakamura
Original assignee: The Ritsumeikan Trust; Toyota Jidosha Kabushiki Kaisha
Priority date: 2004-08-25
Filing date: 2005-08-25
Publication date: 2006-03-02
Also published as: JP4146877B2; US7586440B2; EP1793238A1; DE602005014371D1; EP1793238A4; WO2006022318A8; JPWO2006022318A1; CA2578018A1; EP1793238B1; US20080258966A1; KR101151782B1; CA2578018C; KR20070059109A; CN101014874A; CN101014874B

Abstract

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

Description

明細書

単独測位装置および単独測位方法

技術分野

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

背景技術

[0002] 従来、測位衛星力送信される測位用信号を用いて単独測位を行う装置が各種開示されている。これらの装置による基本測位演算は、コード擬似距離を用いて、受信機の 3次元の位置 (誤差)と受信機の時計誤差とを未知数とした非線形連立方程式に-ユートン法や拡張カルマンフィルタを適用することで求めるものであった。この際

、電離層遅延による誤差や対流圏遅延による誤差の影響を除去する方法が各種考案されている。その 1つとして、これらの遅延による影響を初期状態から「0」に設定する方法があり、またその 1つとして、 1重位相差等の演算処理を行ってこれらの遅延量を推定する方法があった。

非特許文献 1 :日本測地学会編著，「新訂版 GPS—人工衛星による精密測位システム―」，社団法人日本測量協会， 1989年 11月 15日， p. 121 - 140

発明の開示

発明が解決しょうとする課題

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

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

課題を解決するための手段

[0005] この発明は、複数の測位衛星力送信される測位用信号を用いて得られる複数の測位衛星のそれぞれと受信機との距離力ゝら該受信機の位置を測位する単独測位装 ¾【こ; i l /、て、

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

[0006] また、この発明は、複数の測位衛星力送信される測位用信号を用いて得られる複数の測位衛星のそれぞれと受信機との距離力ゝら該受信機の位置を測位する単独測位方法であって、

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

該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算することを特徴として、る。

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

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

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

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

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

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

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

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

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

発明の効果

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

[0017] また、この発明によれば、 LAMBDA法等の整数値バイアス推定法を用いることで

、より高精度に単独測位を行うことができる。

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

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

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

図面の簡単な説明

[0021] [図 1]本発明の実施形態の単独測位装置の概略構成を示すブロック図である。

[図 2]本発明の実施形態の単独測位装置と GPS受信機とからなる単独測位システムの処理フローを示すフローチャートである。

[図 3]本発明の実施形態の単独測位方法を用いた場合の受信機位置 (経度、緯度）を示した図である。

[図 4]本発明の実施形態の単独測位方法を用いた場合の受信機位置 (楕円体高)を示した図である。

符号の説明

[0022] 10—単独測位装置

11—航法メッセージ解析部

12—衛星情報処理部

13—測位演算部

20— GPSアンテナ

30— GPS受信機

発明を実施するための最良の形態

[0023] 本発明の実施形態に係る単独測位装置について図を参照して説明する。なお、以下の説明では、 GPSについて説明する力他の全ての GNSS (全地球的航法衛星システム）に適用することができる。

図 1は本実施形態の単独測位装置の概略構成を示すブロック図である。また、図 2は本実施形態の単独測位装置と GPS受信機とからなる単独測位システムの処理フローを示すフローチャートである。

[0024] 図 1に示すように、単独測位装置 10は、 GPS受信機 30に接続し、航法メッセージ解析部 11、衛星情報処理部 12、および測位演算部 13を備える。

[0025] GPS受信機 30はアンテナ 20に接続し、アンテナ 20で受信した GPS衛星（測位衛星)からの GPS信号 (測位用信号)より、既知の方法で L1波および L2波のキャリア位相、 CZAコードおよび P (Y)コードのコード擬似距離 (擬似距離)を取得するともに、 L1波に重畳された航法メッセージを取得する（sl→s2)。そして、 GPS受信機 30は航法メッセージを航法メッセージ解析部 11に出力する。また、 GPS受信機 30はキヤリア位相およびコード擬似距離を測位演算部 13に出力する。

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

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

[0028] 測位演算部 13は、入力された各観測値、すなわち、選定された各 GPS衛星に対する L1波のキャリア位相、 L2波のキャリア位相、 CZAコードによるコード擬似距離、 Ρ (γ)コードのコード擬似距離、各 GPS衛星の衛星軌道情報、各 GPS衛星の衛星時計誤差情報、電離層遅延量情報、対流圏遅延量情報を用いて、後述する線形回帰方程式を形成する。そして、測位演算部 13は、この線形回帰方程式にパラメータ推定アルゴリズムである最小二乗法を適用させて、 Ll、 L2波に対する整数値バイァス N , N と受信機の位置とを推定演算する（s4)。この演算は推定値の変化が予め

LI L2

設定された所定の閾値以下になるまで繰り返し行われ、推定値の変化が所定閾値に達した時点での受信機位置の推定演算結果が測位結果として出力される。

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

[0030] 次に、前述の整数値バイアス N , N および受信機位置 uの推定演算アルゴリズム

LI L2

について詳述する。

[0031] 一般に、受信機 u、 GPS衛星 pに対するキャリア位相 φ ^ρ の観測方程式は式（1)でし， u

表すことができ、コード擬似距離 (擬似距離） ^Ρ の観測方程式は式 (2)で表される

c，u

。ここで、マルチパス誤差は微少として無視する。

[0032] [数 1]

(0 = { (t, t - ) - 0lS (t) + δ (t) + c[bt _u (t) - 6t^p (t - )J+ N + ε (t) — ( 1)

Pc^P (0 = (t, t ） + 5I^ (t) + δΤ (t) + c[6t_u (t) - 5t^p(t - )J+ e (t) — (2) [0033] ここで、 λ」ま L波の波長を示す。また、 (t, t— τ は時刻 tでの受信機 uと時刻（ t— τ^Ρ )での GPS衛星との距離を示す。 δ Ι^Ρ (t)は Ll波の電離層遅延量を示し、 δ T^P (t)は Ll波、 L2波の対流圏量遅延を示す。 S t (t)は真の時刻 tでの受信機 uの時計誤差を示し、 δ t^P (t— τ ^P )は時刻 (t τ ^P )での GPS衛星 pの時計誤差を示す。 N^P は L波における受信機 uと GPS衛星 pとの間の整数値バイアスを示し、 ε ^Ρ (t) し，し，

, e^P (t)はそれぞれ観測雑音を示す。

c，u

[0034] したがって、 L1波におけるキャリア位相 φ^Ρ の観測方程式は式（3)となり、 L2波に

Ll,u

おけるキャリア位相 Φ ^P の観測方程式は式 (4)となる。

[0035] [数 2]

Φ _u (t) = r^- {r_u ^p (t,t- )- 0I^(t) + δΤ (t) + c[0t_u(t) - St - (t - τζ )}+ N _lu + (t) — (3)

φ _u (t) = ^ fr_u ^p (t， t ） - (t) + δΤ (t) + c[6t _u (t) -6t^p(t-

+ _{2iU +}e^_u(t)

[0036] ここで、 f は Ll波の周波数、 f は L2の周波数を示す。

Ll L2

[0037] さらに、 Φ^Ρ , Φ^Ρ をそれぞれ次式 (3'), (4')により定義する。

Ll,u L2,u

[0038] [数 3]

= (t, t - ） - 6I_U ^P (t) + 0T (t) + c[0t_u (t) -0t"(t- )]+ λυΝ^ + λ_ίιε (t)

Φ λ_ί2φ (t)

= r_u ^p(t,t- τ - (t) ₊ ST (t) + c (t)-St^p(t- τ + _L2N _{U + L} , _u(t)

[0039] また、 CZAコードによるコード擬似距離 ^p の観測方程式は式 (5)となり、 Pコ、

CA

ドによるコード擬似距離 ^P の観測方程式は式 (6)となる。

P，u

[0040] 画

P t) = (t，t -て ₊ &I t) ₊ ST (t) + 4&t_u(t)- St^p(t— " j+e , _u(t) "(5)

P (t) = (t, t - ) + (t) + δ (t) + c[bt (t) - St (t - - )] + e (t) - (6) [0041] ところで、受信機と GPS衛星との距離 r^P (t, t— τ ^P )は、式 (7)で表すことができる

[0042] [数 5]

r (t)_≡r_u ^p(t,t- ) __(?) = V( , (t) - X ^p (t - ))² + (y _u (t) - y " (t - ))² + (z_u (t) - z ^p (t - x_u ^p ))²

[0043] 次に、未知数である受信機位置 u(t)≡[x (t), y (t), z (t)]^Tを先験的な推定受信機位置ョ^ ， y^t), Z_u ^G)(t)]^Tのまわりで 1次のティラー級数展開を行い、 r^P (t)を線形近似すると次式を得られる。

[0044] [数 6]

P = 1， 2 , n

ここで、

',ω 6r

g (9) p = l , 2， · · · , n _s

[0045] この式（9)により、式（3')、 (4，）、 (5)、（6)は、次式（10)〜（13)で表される。

[0046] [数 7]

¾ί = + g^p'⁽ⁱ⁾(u一 υ^ω) - bll + δΤ + c(6t_u -6t^p) + _L1NL_>U + _L1 - (10)

= )(u- u^(j)) δ δΤ + -6t^p) _{+ L2}NP_¾u +λ_ί2ε "(11)

p^)_u = r^p _(i) +g^p'⁽⁾⁾(u-u⁽ⁱ⁾) + 6Ij; +6T_U ^P +c(6t_u -bt^p) + e_C ^p _Au —(12)

. _gP,(j)_(u __u(j)_{) +} liL_{SI +δ +c(0tu}— st^p) + e -(13)

[0047] ：で、 0>_e ^p'w , , e^P,u 、を式 (14)〜（17)で定義する。

[0048] [数 8] Φε^' ≡ -(r^p _(i) -g^p'⁽ⁱ⁾u^(j)) - (14)

Φ_ε ¾ _≡ - (r - '")u('.)) - (15) n_pP.(j)一 _nP,(j) _rP , „P.(j),,(j) —

P^eCA,u = PCA,„― ^r _u(i> + g ^U 、丄 ^W

[0049] これにより、式（10)〜（13)は次式（18)〜（21)となる,

[0050] [数 9]

φ gp'(j)_u _ δΐ^ + δΤ + c(6t_u -0t^p) + k_uN^p _{U u} + _{U u} (18) f

Φε¾^ ≡ g^w"u— ~ +6T_u ^p +c(5t_u -^) + _L2N _2u + _L2e[ (19)

g^p'^(j)u + 0lP +5T_U ^P +c[bt -6t^pj+eP_{A u} —(20)

+_6XuP _+c[st_u -&t^p]+e -(21)

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

[0052] ここで、を式（22)で定義する。

[0053] [数 10]

ι,ϋ·)

2,(j)

G ( _ (22)

[0054] なお、は式（23)で表される

[0055] [数 11] δ > δ

G (23)

& ^δ

δτ^, ( 4^η→· (y^ -y^p) ， (4^j) -_Z ^p、

(ノ） δζ (j)

( P = 1， 2 η _s)

[0056] また、 Φ ®

Ll 、を行列式で表す。

[0057] [数 12]

(24)

[0058] この式を用いることで、式（18)〜式（21)はベクトル行列からなる次式で表される [0059] [数 13]

Φ '> P _u G?

[0060] もしくは、

[0061] [数 14] = H⁽ⁱ⁾e_u + V (26)

c6t_u c6t^p 6I_U 0T_{U u} f _ (29)

[0062] ところで、前述のように、航法メッセージには GPS衛星毎の時計誤差情報が含まれており、航法メッセージ解析部 11で GPS衛星毎の時計誤差情報が取得されて、る。この衛星時計誤差の観測値 δ te^Pは次式で表され、前述の回帰方程式に加えることができる。

[0063] [数 15]

6te^p = [0 0 1,— 0 0 0 0je_u + e_fttP

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

[0065] [数 16]

61e_u = [θ 0 0 I„ 0 0 ojp_u + e_61u

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

[0067] すなわち、式（26)に示した回帰方程式は、次式に置き換えられる。

[0068] [数 17]

Η '°Θ., +υ„ 一（31)

y⁽i≡ [ 4 Φ4¾„ Ped' u Ste^p 01e_u 6Te_u] 一（32)

1 - I - I I I 0

_G ， f²

1 -I - - I I 0 I

f_L2

G(_U ^J) 1 - I I I 0 0

H_E ^(j>≡ f² 一（33)

1 -I I 0 0

f_L ² ₂

0 0 I 0 0 0 0

0 0 0 I 0 0 0

0 0 0 0 I 0 0

υ",Ε ^≡ [ ½ ^L2^£L2,u ^eCA,u ^eP,u ^eSf ^eSI_a ^e6Tu \ 一（34)

[0069] ここで、 υ の共分散行列を Rは式（35)で表される。

uE

[0070] [数 18]

R≡ Cov[v_{u E} ] = (35)

ここで、

≡^λΐ^συ' ≡

[0071] 従って、 θ "の推定値 Θ e⁰^は、式（36)となる。 [0072] [数 19]

6e^(j+1) = {Η^') ]^— ' [H^ ]}— '[ f R- - (36)

[0073] また、推定値 Θ e⁰^の分散値は、式 (37)となる。

[0074] [数 20]

Cov[e_e叫ョ E[(ee⁽ⁱ⁺¹⁾ - 9)(6e^(j'⁺¹⁾ - θ)^τ

- (37)

[0075] この式（36)に基づ、て、最小二乗法を用いて繰り返し推定演算を行うことで、 Θを構成する未知数を収束させて各値を取得する。この際、説明変数ベクトル Θを構成する未知数の収束条件は、未知数の推定値の繰り返し演算での変化量が、予め設定した閾値以下に達した時点である。そして、この条件下で説明変数ベクトル Θの構成要素の値を確定させる。このような演算を行うことで、整数値バイアス N , N

LI L2ゃ受信機位置 uを推定することができる。

[0076] このように、本実施形態の構成を用いることで、 1重位相差等の処理を行うことなく整数値バイアス N , N や受信機位置 uの推定演算、すなわち単独測位を行うことが

LI L2

できる。これにより、従来例より簡素なアルゴリズムで単独測位を行うことができ、且つ単独測位の演算処理速度を向上することが可能になる。

[0077] ところで、このような推定演算を用いる場合、未知パラメータは Θ の各要素であり、式 (31)に示す回帰方程式内の未知パラメータ数は、受信機位置の uが 3個、受信機時計誤差 δ t力 Si個、 GPS衛星の時計誤差 δ t^Pが η個（GPS衛星数)、電離層遅延

u a

δ Iが η個、対流圏遅延 δ Τが η個、 L1波整数値バイアス Ν 力 ¾個、 L2波整数値 s s し 1 s ノィァス Ν 力個、の合計 4 + 5η個となる。一方で、観測値は、 L1波のキャリア位

L2 s s

相、L2波のキャリア位相、CAコード擬似距離 p e® 、Pコード擬似

Ll,u L2,u CA,u 距離 p 、 GPS衛星の時計誤差 δ te 電離層遅延 δ Ie、対流圏遅延 δ Teがそ

CA

れぞれ測位可能な GPS衛星数に相当する nであり、合計 7n個となる。

[0078] 従って、この回帰方程式で前記未知数の解を確定する成立条件は、次式となる。

[0079] 4 + 5n≤7n

s s

これは n≥ 2を意味する。すなわち、少なくとも 2つの GPS衛星力もの測位用信号を s

受信すれば受信機の単独測位が可能となる。し力しながら、この回帰方程式を最小二乗法で解く場合、式 (33)でのの逆行列が必要となる。ところが、測位衛星数 n 力^の場合には前記行列が特異行列となり、 nが 3の場合には前記行列が特異となる可能性が非常に高くなり、逆行列が存在せず現実的に演算できない。

[0080] また、前述の実施形態では、 Pコードを用いて推定演算を行ったが、 P (Y)コードは秘匿コードであるため、コード擬似距離の観測は事実上難しい。この場合、観測値が 1種類 (η個)減少するので回帰方程式の成立条件は次式となる。

[0081] 4 + 5η≤6η

これは η≥4を意味し、実質的に少なくとも 4つの測位衛星力の測位用信号を受信すれば受信機の単独測位が可能となる。すなわち、少なくとも 4つの測位衛星からの測位用信号を受信できれば、前述の単独測位は可能となる。

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

[0083] 図 3は本実施形態の単独測位方法を用いた場合の受信機位置 (経度、緯度)を示した図である。図 4は本実施形態の単独測位方法を用いた場合の受信機位置 (楕円体高）を示した図である。

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

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

[0086] なお、このシミュレーションでは、受信機位置の初期値として、 RINEXデータの" A PPROX POSITION XYZ"に記載された座標を用いた。また、電離層遅延量は、航法メッセージに含まれる電離層遅延に関するデータから、いわゆる放送モデル (K1 obucharモデル)を適用して算出した値を用い、対流圏遅延量は、航法データから計算される測位衛星の仰角情報力次式を適用して算出した値を用いている。

[0087] [数 21]

6Te = ²'⁴⁷ [ m ]

sin ξ + 0.0121

[0088] ここで、 ξは仰角である。

さらに、観測誤差の分散は、それぞれ、コード擬似距離が 1. 5 [m] (1 σ )で、キヤリァ位相が (波長 λ ΖΐΟ) + (1. 5/10) [πι] (1 σ )で、衛星時計誤差が 3. 6 [m] (l σ )で、電離層遅延が 7. 0[m] (l σ )で、対流圏遅延が 0. 7[m] (l σ )でそれぞれ設定されている。また、このシミュレーションは、最小二乗法による繰り返し演算での、受信機位置の推定値の変化のノルム値が 1 X 10— ³[m]以下になるまで計算を繰り返した。

[0089] [表 1]

[0090] [表 2]

[0091] [表 3] 緯度経度楕円体局

平均値 34. 98186963 [度] 135. 9641 1 15 [度] 236. 1081403 [度] 標準偏差 1 . 347937297 [m] 1. 572959303 Cm] 3. 184914299 M 相対測位結果 34. 981876 [度] 135. 964098 [度]

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

[0093] 以上のように、本実施形態の単独測位装置および単独測位方法を用いることにより

、従来よりも簡素な演算処理で高精度に受信機の単独測位を行うことができる。

[0094] なお、前述の最小二乗法による推定演算時に、式 (37)より得られる分散値を用いて、 LAMBDA法を適用することで、整数値バイアス N , N をフィックスして、さらに

LI L2

高精度の受信機位置を測位することができる。

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

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

[0097] 6n XK=6Kn

s s

一方、観測エポックが複数となるキネマティック測位では、整数値バイアス N , N

LI L2 を除く各未知数 (未知パラメータ）がエポック数に応じて増加するので、未知数の数は、次式で表される。

[0098] 4 + 5n + (K— 1) (4 + 3η ) =4Κ+ (3Κ+ 2) η

s s s

したがって、回帰方程式を最小二乗法により解くことができる必要条件は、次式となる。

[0099] 6Kn≥4Κ+ (3Κ+ 2) η

s s

(3K- 2) n≥4K すなわち、

n≥4K/(3K-2)

そして、エポック数 Kを 2→∞の場合、 nの取り得る整数値を計算すると、 GPS衛星数 nは次式の値となる。

[0100] n≥2

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

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

[0102] まず、式 (31)より、エポック数に対応する観測時間 tを要素として表現し直すと、 [0103] [数 22]

yl!⁾ _E(t) = H_E ^(j)(t)6_u(t) _{+ Uu>E}(t) -(38)

t = 1 , 2， 3， · · ·

[0104] これを、ベクトル行列表示すると、式（39)になる。

[0105] [数 23]

1 一 I - 1

Φ (り f ²

1 -I

f²

ρ ) 1 -I 1

p '_u(t) (39)

1 一 I I

c6te^p(t)

0 0 I 0 0 0 0

5Ie_u(t)

0 0 0 I 0 0 0

6Te_u (t)

0 0 0 0 1 0 0

t = 1 , 2 3 , ■

[0106] ここで c δ t^P, δ I , δ Tをそれぞれ式（40a) , (40b) , (40c)とする。

[0107] [数 24]

c5t^p(t) = c6te^p(t)-e_ft„(t) —(40a)

6I_u(t) = 5Ie_u(t)-e₆₁ (t) _ (40b)

6T_u(t) = 6Te_u(t)-e_6T (t) —(40c) [0108] これらを式（39)の回帰方程式に代入し、さらに、次式に示す定義を行う _c

[0109] [数 25]

y¾_u (t) = e^_u (t) + c6te^p(t) + 6Ie_u (t) - 5Te_u (t) - (41)

y (t) = e¾,_u (t) + c0te^p (t) + -^6Ie_u (t) - 0Te_u (t) - (42)

(t) = pe _u (t) + c0te^p(t)― 01e_u (t) - 0Te_u (t) - (43)

(t) = pe (t) + c6te^p(t) - - 6Ie_u (t) - STe_u (t) _ (44)

[0110] この結果、新たな回帰方程式は、式 (45)で表される。

[0111] [数 26]

u

(45)

[0112] ここで、未知の受信機位置 u (t)の速度に対するマルコフ過程モデルを仮定し、かつ、受信機の時計誤差 c δ t (t)に対してマルコフ過程モデルを仮定できるものとし、速度ベクトルを V (t)とおき、受信機の新たな状態ベクトル 7? (t)を次式で定義する。

[0113] [数 27]

u(t)

v(t)

"_u (t) c6t_u (46)

N" ,,

[0114] そして、この状態ベクトル r? _u (t)について次の状態方程式を構成する。

[0115] [数 28]

Tl_u (t + 1) = Ari_u (t) + Bw(t) —（47)

[0116] この際、式 (45)の観測方程式は、次式で表すことができる。 [0117] [数 29]

一（48)

[0118] この観測方程式をベクトル表示すると、式 (49)となる。

[0119] [数 30]

yl⁾ _R (t)

_{+ UuiR} (t) - (49)

t = 1 , 2， . . . ·

[0120] これは、すなわち、式 (47)、式 (48) (式 (49) )力もなるカルマンフィルタを意味する。この際、観測誤差 υ (t)の共分散行列 Rは、式 (35)と同様に求めることができる。そして、これら状態方程式 (47)、観測方程式 (48)、を用いることにより移動する受信機位置を推定演算することができる。

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

[0122] ところで、前述の説明では、 GPS衛星の位置を推定演算しな力つた力次に示す方法 (アルゴリズム）を用いることにより、 GPS衛星の位置を推定演算することができる

[0123] 式 (7)に示すように、受信機と GPS衛星との距離を定義し、受信機とともに GPS衛星の位置を推定演算する場合、受信機位置に関する線形近似式は、前述のように表される。

[0124] [数 31]

(z_u - z^p )

( p = 1 , 2， · · · ， n _s )

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

[0126] [数 32] _ (z„ - Z^p)

r_u ^p

( ρ = 1， 2 , · · · , η _s )

[0127] そして、 sp≡[x^P, y^P, z^P]^Tと定義すると次式の関係が成り立つ。

[0128] [数 33]

6r_u ^p _ _ 0r_u ^p

0u 6s^p

[0129] したがって、 u, spの推定値は、それぞれ ι se^Pと 1次ティラー級数展開とを用いて次式で近似される。

[0130] [数 34]

[u-s^p _(u ('に se^p)] -(50)

διι

( P = l , 2， · · ·， n _s)

[0131] したがって、式（10)力式（13)は、それぞれ次式で表される。

[0132] [数 35]

Φ = +g^p'⁰⁾[u-s^p-(u^<i) -se^p)J-6IJ ₊δΤ_υ ^ρ ₊ο(δΐ_υ -δΐ") ₊ λ,₁Ν^ ₊λ_ί1ε _1ιϋ 一（51)

(52) + c(6t_u _&t^p) + e

_A (53)

■ gP'")[_u _: s^p -(u⁽ⁱ⁾ - se^p)J+ +ST_u ^p +c(St_u -6t^p) + _eP (54)

[0133] ここで、 , e^P,u 、をそれぞれ式 (55)〜（58)で定義

する。

[0134] [数 36] Φ _≡φ に [ -g^p'⁽ⁱ⁾(u⁽ⁱ⁾ -se^p)J -(55)

*e^ ^^ -fc» -g^P'⁽ⁱ⁾(^u(i) -^seP)J 一（56)

p ^_u≡ P ― [r_u ^p ₍₁₎ - g^p'⁽ⁱ⁾(u⁽ⁱ⁾ - se^p)J — (57)

pe^p _f ≡ p^p _P,(_u ^j> - [r_u ^p _(i) - g^p'⁽ⁱ⁾(u ^j) - se^p)J 一（58)

[0135] これにより、式（51)〜（54)は次式（59)〜（62)となる。

[0136] [数 37]

Φε ^≡ g^p'⁽ⁱ⁾u― g^p'('V― 61；； + 6T_U ^P + c(6t_u -6t^p) + _L1N _u + Λι,η ― (⁵⁹)

_Φε 】½— -1 L_6IP _{+δτ +C}(_6tu - δ + _L2N + _{L2e 2},_u -(60) e^^D _u +0T_U ^P +c[6t_u -&t^p_|+ _Au —(61) f ² _c r\ Ί

gp.(j)_u .U J_5]P _{+ 6T}P ₊ —_5tpj_{+ e} -(62)

上 2

[0137] ここで、 Φ ® , Φθ^Ρ,(ί) , e^P,(i) , θ^Ρ,α) 、をそれぞれ次式で定義する。

Ll,u L2,u Ll,u L2,u

[0138] [数 38]

[0139] これにより、式（59)〜式（62)はベクトル行列からなる次式で表される。

[0140] [数 39]

ただし、行歹 IJ )

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

[0142] [数 40]

sk = s + e _s，

，A，ひ¹)⁷"]⁷"， 3n_s x l

[0143] これを式 (63)に適用すると式（64)となる。

[0144] [数 41]

一（64)

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

[0146] さらに、以下のアルゴズムを用いることで、このような場合でもカルマンフィルタを適用することができる。 [0147] すなわち、 s， c δ t^p, δ I_u, δ をそれぞれ式 (65a)〜（65d)で定義する, [0148] [数 42]

s(t) = sk(t)-e_s (f) 一（65a) c6t^p(t) = c6te^p(t)-e₈,_p (t) —(65b)

5I_u(t) = 5Ie_u(t)-e_6Iu (t) —(65c)

6T_u(t) = 0Te_u(t)-e_6T (t) — (65d)

[0149] これらの関係を用い、式 (64)の回帰方程式に代入し、さらに

次に示す定義を行う。

[0150] [数 43]

y (0 = Φ4ί)" ( + G (t)sk(t) + cdteP (t) + &e_u (t) - STe_u (t) - (66) g„ ） (t) - STe_u (t) (67)

y¾_u (0 =

(0 + G⁽ (t)sk ( + cdteP (t) - 6Je_u (t) 6Te_u (t) (68) ¾ ( = pe_P ^(j (0 + G¾ (t)sk(t) + cdteP (t) - 6Ie_u (t) - 5Te_u (t) ― (69)

[0151] この結果、新たな回帰方程式は、式（70)で表される。

[0152] [数 44]

y _u (t)

(t)

G⁽»_ue_s (t) + e_fttP (t) + c (t)一 e_6Tu (t) + _L1e_L1>u

G ) (t) 1 I 0 u(t)

G« e_s(t) + e_6tP (t) + ^ Le_¾ (t)一 e_8To (t) + _L2E_L

G⁽ _u ^j) (t) 1 0 I

G<^j) (t) 1 0 0

G« e_s (t) + e_S(P (t)一 e_61u (t)一 e_6Tu (t) + e_{CA u}

Gi^J) (t) 1 0 0 λ,，： N,つ

G^_ue_s (t) + e_6(P (t)― e_6Iu (t) - e _w' (t) + e_{P u}

" (70)

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

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

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

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

[0157] また、前述の実施形態では、 LAMBDA法を利用する例を示した力他の整数値バイアス推定方法を用いてもょ、。

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

Claims

請求の範囲

[1] 複数の測位衛星力送信される測位用信号を用いて得られる前記複数の測位衛星のそれぞれと受信機との距離力該受信機の位置を測位する単独測位装置において、

前記測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値カゝら前記測位衛星の軌道情報および衛星時計誤差を観測する衛星情報検出手段と、

電離層遅延量情報を取得する電離層遅延情報取得手段と、

対流圏遅延量情報を取得する対流圏遅延情報取得手段と、

受信機位置を過去の受信機位置推定結果および前記測位衛星の軌道情報を用いて線形近似化し、

未知数である線形近似された受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とする回帰方程式を構成し、該回帰方程式にパラメータ推定ァルゴリズムを適用することで、少なくとも受信機位置を推定演算する測位演算手段と、を備えたことを特徴とする単独測位装置。

[2] 前記測位演算手段は、整数値バイアス推定法を用いて前記整数値バイアスをフィッタスして、前記受信機位置を推定演算する請求項 1に記載の単独測位装置。

[3] 複数エポックで取得した前記キャリア位相、前記コード擬似距離、前記受信機時計誤差、前記衛星軌道情報、前記衛星時計誤差、前記電離層遅延量情報、および、前記対流圏遅延量情報を記憶する一時記憶手段を備え、

前記測位演算手段は、取得した前記複数エポックの各情報を用いて、前記回帰方程式を構成する請求項 1または請求項 2に記載の単独測位装置。

[4] 前記測位演算手段は、前記回帰方程式に状態推定アルゴリズムを適用する請求項 1〜3のいずれかに記載の単独測位装置。

[5] 前記測位演算手段は、前記衛星軌道情報を目的変数に含み、前記衛星軌道情報に対応する衛星軌道誤差を説明変数に含む前記回帰方程式を用いて、測位衛星位置を推定演算する請求項 1〜4のいずれかに記載の単独測位装置。

[6] 複数の測位衛星力送信される測位用信号を用いて得られる前記複数の測位衛星のそれぞれと受信機との距離力該受信機の位置を測位する単独測位方法であつて、

前記測位用信号に含まれる航法メッセージまたはオフライン処理により予め推定された値カゝら前記測位衛星の軌道情報および衛星時計誤差を観測し、

電離層遅延量情報、および対流圏遅延量情報を取得する対流圏遅延情報取得し過去の受信機位置推定結果および前記測位衛星の軌道情報を用いて受信機位置を線形近似化し、

未知数である線形近似された受信機位置、整数値バイアス、受信機時計誤差、衛星時計誤差、電離層遅延量、および対流圏遅延量を説明変数とし、観測値であるキャリア位相、コード擬似距離、衛星時計誤差、電離層遅延量情報、および対流圏遅延量情報を目的変数とする回帰方程式を構成し、

該回帰方程式にパラメータ推定アルゴリズムを適用することで、少なくとも受信機位置を推定演算する、

ことを特徴とする単独測位方法。

[7] 整数値バイアス推定法を用いて前記整数値バイアスをフィックスして、前記受信機位置を推定演算する請求項 6に記載の単独測位方法。

[8] 複数エポックで取得した前記キャリア位相、前記コード擬似距離、前記受信機時計誤差、前記衛星軌道情報、前記衛星時計誤差、前記電離層遅延量情報、および、前記対流圏遅延量情報を一時記憶し、

取得した前記複数エポックの各情報を用いて、前記回帰方程式を構成する請求項 6または請求項 7に記載の単独測位方法。

[9] 前記回帰方程式に状態推定アルゴリズムを適用する請求項 6〜8の、ずれかに記載の単独測位方法。

[10] 前記衛星軌道情報を目的変数に含み、前記衛星軌道情報に対応する衛星軌道誤差を説明変数に含む前記回帰方程式を用いて、測位衛星位置を推定演算する請求項 6〜9の、ずれかに記載の単独測位方法。