第1章

RTK-GPS の原理と応用
Principle and Practice of RTK-GPS

浪江 宏宗（防衛大学校）

1.1 はじめに本章では、RTK-GPS の測位原理について、筆者が実際に取得したデータを交えて詳述する。まず、

図 1.1 に RTK-GPS の構成概念図を示す。RTK-GPS の主な構成要素は、既知位置に設置した固定基準局（ Reference Station, Base Station 以下基準局と記述する ）、データ伝送システム、および利用者局（ユーザ局： User Station, 移動局： Remote Station, Rover Station ）である。既知の固定点に設置された基準局の受信機では、常時可視 GPS 衛星からの電波を受信し、搬送波位相の積算値データ（ Carrier

Phase 以下、搬送波位相データ）を測定する。そして、それらのデータを含む測位用データを利用者局に伝送する。利用者局でも同様に搬送波位相を測定し、基準局から伝送されてきた測位用データと共に使用して、利用者局の三次元位置が求められる。

図 1.1: RTK-GPS の構成概念図

GPS 受信機で取得したデータは、衛星からの搬送波の、波数を単位とした位相データである。そこで、
まず搬送波位相データについて述べる。次に、mm オーダの高精度測距を実現する、搬送波の一重位相差および二重位相差について記述する。そして、搬送波位相データのアンビギュイティ（ Ambiguity1 ）、
および RTK-GPS の測位計算について述べる。さらに、最後に RTK-GPS 測位の幾何学的な解釈について記述する。
1 整数値バイアス。搬送波位相データを測距に使用したときに、必然的に発生する搬送波長の整数倍の位相不明確量。

1.2 搬送波位相データおよび位相差(1) 搬送波位相の積算値データ
RTK-GPS では衛星からの搬送波位相積算値データを用いて測位計算を行う。電離層や対流圏等による誤差、および受信系の S/N 比等による搬送波位相データの観測の偶然誤差を無視すると、受信機 A で測定された衛星 i からの搬送波位相データ iA は、式 (1.1)(1) で表される。
fiA = riA + di + A + NiA
c

(1.1)

ただし、

f
c
riA
di
A
N

：
：

搬送波周波数光速：

GPS 受信機アンテナ A と衛星 i 間の真の距離
衛星 i の時計誤差による搬送波位相データの誤差
受信機 A の内部時計誤差による搬送波位相データの誤差：

アンビギュイティ（整数値）

：
：

である。

図 1.2 に実際に NovAtel Inc. 社製 GPS 受信機 RT-20（以下 NovAtel RT-20 ）から得た、搬送波

位相データの一例を時系列で示す。同図には複数の衛星から同時に得たデータを一括して示した。
搬送波位相データの増減は緩やかであるように見える。しかし、縦軸は搬送波の波数を単位として搬送波位相データ（波数）×10

7

107 のオーダであり、大きく増減していることが分かる。

–11

–12

–13
0

10

20
30
40
時 間（分）

5060

図 1.2: 搬送波位相データの時系列（複数衛星）
次に、図 1.3 に SA による搬送波位相積算値データの、揺らぎ誤差の一例を時系列で示す。同図から、平均 60 m 程度の振幅で変動していることが分かる。同データは、2000 年 5 月の SA 解除前に取得したものであるため、SA の影響が出ている。SA は衛星内の基準発振器の位相に対して施されていることが分かる。SA 解除後のデータでは、まだ解析していない。

(2) 衛星間一重位相差

次に、ある時刻に GPS 受信機 A で受信している衛星群中の 2 機を、それぞれ衛星 i 、および j とする。また、測定された搬送波位相データをそれぞれ iA , jA とすると、受信機 A に対する搬送2

搬送波位相の揺らぎ（m）

80
40
0
–40
–80
0

10

20
時

30
40
間（分）

50

60

図 1.3: 搬送波位相データの揺らぎ誤差の時系列
波位相データの衛星間一重位相差 DijA は、

DijA = jA 0 iA(1.2)

で表される。2 つの搬送波位相データ iA , jA に共通して含まれていた、受信機の時計の誤差による搬送波位相誤差が、この処理によって完全に消去される。しかし、衛星の時計の誤差による位相誤差は残る。

(3) 受信機間一重位相差

ある時刻に、受信機 A, B でそれぞれ測定した衛星 i からの搬送波位相データを iA , iB とすると、衛星 i に対する搬送波位相データの受信機間一重位相差 DiAB は、

DiAB = iB 0 iA(1.3)

で表される。2 つの搬送波位相データ iA , iB に共通して含まれていた、衛星内部の時計誤差による位相誤差が、この計算により完全に消去される。しかし、GPS 受信機内部の時計誤差による位相誤差は残る。

(4) 二重位相差

ある時刻に、衛星 i と j の搬送波位相データを、それぞれ A, B 2 台の GPS 受信機で測定する。

このとき、各受信機で測定された衛星 i からの搬送波位相を iA , iB し、衛星 j からの搬送波位相データを jA , jB とすると、搬送波位相データの二重位相差 DDijAB は式 (1.4)(1) のようになる。

DDijAB = (jB 0 jA ) 0 (iB 0 iA )= (jB 0 iB ) 0 (jA 0 iA )(1.4)

ここで、GPS 受信機内部の時計誤差、および衛星内部の時計の誤差による位相誤差、SA による

位相誤差、さらに電離層や対流圏に起因する位相誤差が消去される。衛星 i, j と点 A の位置が既知なので、DDijAB を求めることができれば jB

0  、すなわち未知の点 B から衛星 i, j までのiB

距離差に関連した値2 が測定できる。以上のように、RTK-GPS では二重位相差の処理を行うこと
で、SA 等による搬送波位相誤差を消去して、mm オーダの高精度測距を実現している。

2 整数値アンビギュイティを決定できれば距離差と同義

3図 1.4 に前述の搬送波位相データを使用して計算した、二重位相差の一例を時系列で示す。同図のように二重位相差はサイクル スリップ 3 が発生していなければ、1 波長以上の位相のジャンプは二重位相差（波数）なく、連続的に変化する。

1
0
–1
–2
–3
0

10

20
30
40
時 間（分）

5060

図 1.4: 二重位相差の時系列

測 距 誤 差（mm）

80
40
0
–40
標準偏差= 8.8 mm
–80
0

10

20
時

30
40
間（分）

50

60

図 1.5: 二重位相差の測距誤差時系列
さらに、図 1.5 にこの二重位相差の回帰曲線からの残差、すなわち、ランダムな測距誤差の一例を時系列で示す。同図より、二重位相差の計算によって相殺されなかったランダムな誤差は、回帰曲線に対して標準偏差が 8.8 mm 程度であることが分かる。これが NovAtel RT-20 を使用した

場合の、搬送波位相データによる測距精度であるといえる。この値に DOP4 値を乗じたものが、

RTK-GPS の測位精度の性能限界となる。

3
4

Cycle slip： 衛星からの電波が遮断されることにより、アンビギュイティの値が変わること。
Dilution of Precision： 衛星配置にのみ依存する測位精度劣化指数4

1.3 RTK-GPS 測位計算アルゴリズム - FLOAT 解次に、RTK-GPS の測位計算アルゴリズムについて詳述する。まず、FLOAT 解の計算アルゴリズムについて記述する。

GPS では、中心が地球重心（中心）と一致した楕円体を基準にした、WGS-84（ World Geodetic System
1984： 全世界的測地系 1984 ）と呼ばれる座標系により測位計算が行われる (2) 。そして、私達が日常使用する経度, 緯度、および高さの座標に変換される。また、使用している xyz 三次元直交座標系は、ECEF（ Earth Centered Earth Fixed ）と呼ばれている。これは、地球の中心を原点とし、自転軸方向
を z 軸（天の北極方向）、グリニッジ基準子午面と赤道が交わる方向を x 軸。さらに、これらの 2 軸と右手系をなすように y 軸を設定したものである。
各衛星からの搬送波の周波数が等しい（ド リフト無し ）と仮定すると、前述の二重位相差の式 (1.4)
は、エポック t において次式のようにも書き表すことができる。
ccDDijAB (t) = (rjB (t) 0 rjA (t)) 0 (riB (t) 0 riA (t)) + NijAB (t)ff= (rjB (t) 0 riB (t)) 0 (rjA (t) 0 riA (t)) + c NijAB (t)
f

(1.5)

ここで、r はそれぞれ GPS 衛星と利用者局の GPS アンテナ間の真の距離を表している。

基準局の位置は既知であると仮定して、未知のパラメータを整理する。衛星位置 (xi (t), yi (t), zi (t)),(xj (t), yj (t), zj (t)) は、衛星から送信される航法メッセージによって既知である。したがって、利用者
局の xyz 三次元直交座標成分 (xB (t), yB (t), zB (t)) と、各二重位相差の式毎の整数値アンビギュイティということになる。式 (1.5) の右辺の整数値アンビギュイティは、衛星からの電波を受信し続けていれば変化しない。つまり、二重位相差の各式毎に一定値となり、未知数は各式毎に 1 つずつである。結局、独立な二重位相差から定義される方程式の数よりも未知のパラメータが多いので、このままでは未知数を求めることができない。しかし、別のエポック（時刻）にも同様に異なる方程式が定義されるので、利用者局位置も固定であると仮定すれば、整数値アンビギュイティと共に利用者局位置を求めることができる。

まず、式 (1.4) 、および (1.5) を使用する。衛星 1 を基準衛星5 とし、測位に使用した衛星数 s より 1

つ少ない数、独立な二重位相差を定義する。すなわち、エポック t ((1.6) のように定義される。

8>> c DD12
>> f
>> c
>> f DD13><c> f DD14>>>>>>>: c DD1f 1) において独立な方程式が式(t) + (r2A (t) 0 r1A (t)) = (r2B (t) 0 r1B (t)) + c N12AB (t)f= (r3B (t) 0 r1B (t)) + c N13AB (t)AB (t) + (r3A (t) 0 r1A (t))
f
= (r4B (t) 0 r1B (t)) + c N14AB (t)AB (t) + (r4A (t) 0 r1A (t))
f
AB

sAB

..
.

..
.

(1.6)

..
.

(t) + (rsA (t) 0 r1A (t)) = (rsB (t) 0 r1B (t)) + c N1sAB (t)f

実際の測位計算では、エポック t における計算に、エポック 0 における衛星の位置を使用する。すな

5 独立な搬送波の二重位相差をより容易に計算するために、基準衛星を決めている。

5わち、エポック 0 においても独立な二重位相差の方程式が式 (1.7) のように定義される。

8>> c>> f DD12>> c>> DD13<f
> fc DD14>>>>>>>: c DD1f(0) + (r2A (0) 0 r1A (0)) = (r2B (0) 0 r1B (0)) + c N12AB (t)f= (r3B (0) 0 r1B (0)) + c N13AB (t)AB (0) + (r3A (0) 0 r1A (0))
f
= (r4B (0) 0 r1B (0)) + c N14AB (t)AB (0) + (r4A (0) 0 r1A (0))
f
AB

sAB

..
.

..
.

(1.7)

..
.

(0) + (rsA (0) 0 r1A (0)) = (rsB (0) 0 r1B (0)) + c N1sAB (t)f

ここで、r はそれぞれ GPS 衛星と利用者局の GPS アンテナ間の距離を表しているので、基準局の受信機 A 、および利用者局の受信機 B のアンテナ位置 (xA , yA , zA ), (xB (t), yB (t), zB (t)) 、さらには

衛星 i（ i = 1; 2; 3;

1 1 1 ; s ）のエポック 0 における位置 (x (0), y (0), z (0)) と、エポック t における位置i(xi (t), yi (t), zi (t)) の各 xyz 三次元直交座標を使用すると、ii

8>> r (0) = p(x 0 x (0))2 + (y 0 y (0))2 + (z 0 z (0))2>>>< r (0) = p(x (t) 0 x (0))2 + (y (t) 0 y (0))2 + (z (t) 0 z (0))2>> r (t) = p(x 0 x (t))2 + (y 0 y (t))2 + (z 0 z (t))2>>>:p(x (t) 0 x (t))2 + (y (t) 0 y (t))2 + (z (t) 0 z (t))2
r (t) =
iA

A

iB

B

iA

A

iB

B

i

A

i

i

i

B

A

i

A

i

i

B

i

A

i

B

i

B

i

(1.8)

i

で書き表される。しかし、r はそれぞれ未知数の 2 乗および平方根があって、非線形であるので、線形方程式として解くことはできない。このような場合には、未知数をその近似値と修正値との和で表す。
そして、式をその修正値についてテーラ展開する。ここで、修正値は微小であると仮定して、二次以上の高次項を無視し、式を線形化する。そうすれば、修正値についての連立一次方程式となる。さらに、
最小自乗法による逐次近似計算で、必要な精度まで計算を繰り返し、未知数を求めることができる。

まず、エポック t における n 回目の逐次演算後の利用者局位置 (xB;n (t), yB;n (t), zB;n (t)) とすると、そのときの衛星 1（基準衛星）と利用者局 B 間、および衛星 i と利用者局 B 間の近似距離差 r1iB;n (t) =riB;n (t) 0 r1B;n (t) は、

r1iB;n (t)

=

q
(x
q

B;n

(t) 0 xi (t))2 + (yB;n (t) 0 yi (t))2 + (zB;n (t) 0 zi (t))20 (x (t) 0 x1 (t))2 + (y (t) 0 y1 (t))2 + (z (t) 0 z1 (t))2
B;n

B;n

B;n

(1.9)

で表される。ただし、n = 0 における利用者局の位置 (xB;0 (t), yB;0 (t), zB;0 (t)) は、逐次計算の初めに適当に与えた位置である。また、近似距離差 r1iB;0 (t) は衛星 1 および衛星 i と初めに与えた利用者局 B間の計算距離差である。

さらに、エポック t において測位計算に使用する、エポック 0 の衛星 1 と利用者局 B 間、および衛星
i と利用者局 B 間の近似距離差 r1iB;n (0) = riB;n (0) 0 r1B;n (0) は、

r1iB;n (0)

=

q
(x
q

B;n

(t) 0 xi (0))2 + (yB;n (t) 0 yi (0))2 + (zB;n (t) 0 zi (0))20 (x (t) 0 x1 (0))2 + (y (t) 0 y0 (t))2 + (z (t) 0 z1 (0))2
B;n

B;n

で表される。

6

B;n(1.10)

利用者局の位置、および各二重位相差のアンビギュイティ値 N は初めに適当に与えるので、式 (1.6)と式 (1.7) は、実際は両辺に対して等号は成立しない。そこで、これらの各式において（左辺）-（右辺）B (t) を次のように定める。
0c
BB f DD12 (0) + r12 (0) 0 r12
BB 
BB c DD13 (0) + r13 (0) 0 r13
BB  f


BB c
DD
(0)
+
r
(0)
0
r14
BB f 14
14
BB
..
..
.
.
BB


BB  cB f DD1 (0) + r1 (0) 0 r1
B (t) = BBB BB c DD12 (t) + r12 (t) 0 r12
BB  f


BB c
DD
(
t
)
+
r
(
t
)
0
r13
13
13
BB f
BB  cBB DD14 (t) + r14 (t) 0 r14
BB f
..
..
BB
.
.
B@ 


c

を計算して、行列

n

AB

A

AB

A

AB

A

sAB

sA

AB

A

AB

A

AB

A

..
.

0

DD1sAB (t) + r1sA (t)

CC
C
 CCC
(t) C
 CCC
(t) C
CC
 CC
(t) C
C
 CCC
(t) C
CC
CC
 CCA

(0) + c N1sAB;nfcNB;n (t) +f 12AB;nc
N
B;n (t) +f 13AB;nc
N
B;n (t) +f 14AB;n

sB;n

n

f

1

c
N
(t) C
B;n (0) +CC
f 12AB;n

CC
c
(0)
+
N
(
t
)
CC
B;n
f 13AB;n
 CC
c
N
(t) C
CC
B;n (0) +f 14AB;n..
.

(1.11)cr1sB;n (t) + N1sAB;n (t)f

また、利用者局の近似位置 (xB;n (t), yB;n (t), zB;n (t)) に対する修正値を、それぞれ 1xB;n (t), 1yB;n (t),1zB;n (t) とすると、両者の関係は式 (1.12) のようになる。

8
>> x
<
>> y
:z

(t) = xB;n (t) + 1xB;n (t)= yB;n (t) + 1yB;n (t)B;n+1 (t)= zB;n (t) + 1zB;n (t)B;n+1 (t)B;n+1(1.12)

さらに、アンビギュイティN1iAB;n (t) に対する修正値を 1N1iAB;n (t) とすると、両者の関係は式 (1.13)のようになる。

8N>> 12 +1 (t) = N12>> N>< 13 +1 (t) = N13> N14 . +1 (t) =. N14
>>
..
..
>>
: N1 +1 (t) = N1
AB;n

AB;n

AB;n

sAB;n

(t) + 1N12AB;n (t)AB;n (t) + 1N13AB;n (t)AB;n (t) + 1N14AB;n (t)..
.
sAB;n (t) + 1N1sAB;n (t)AB;n(1.13)

ここで、N1iAB;0 (t) は各エポックにおいて、逐次計算の初めに適当に与えるアンビギュイティ値である。
次に、式 (1.9), (1.10) から、概略の利用者局位置から見た、衛星へ向かう視線の方向余弦の差の、各xyz 座標成分である @r1sB;n (0) , @r1sB;n (0) , @r1sB;n (0) , @r1sB;n (t) , @r1sB;n (t) , @r1sB;n (t) を計算す
@xB;n (0) @yB;n (0) @zB;n (0) @xB;n (t) @yB;n (t) @zB;n (t)7ると次式のようになる。

8
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
<
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
:

(0) = @r1iB;n (0)@xB;n (0)xB;n (t) 0 x1 (0)= p
2(xB;n (t) 0 x1 (0)) + (yB;n (t) 0 y1 (0))2 + (zB;n (t) 0 z1 (0))2
xB;n (t) 0 xi (0)0p2(xB;n (t) 0 xi (0)) + (yB;n (t) 0 yi (0))2 + (zB;n (t) 0 zi (0))2= @r1iB;n (0)i;n (0)@yB;n (0)yB;n (t) 0 y1 (0)= p(xB;n (t) 0 x1 (0))2 + (yB;n (t) 0 y1 (0))2 + (zB;n (t) 0 z1 (0))2
yB;n (t) 0 yi (0)0p2(xB;n (t) 0 xi (0)) + (yB;n (t) 0 yi (0))2 + (zB;n (t) 0 zi (0))2= @r1iB;n (0)i;n (0)@zB;n (0)zB;n (t) 0 z1 (0)= p(xB;n (t) 0 x1 (0))2 + (yB;n (t) 0 y1 (0))2 + (zB;n (t) 0 z1 (0))2zB;n (t) 0 zi (0)0p2(xB;n (t) 0 xi (0)) + (yB;n (t) 0 yi (0))2 + (zB;n (t) 0 zi (0))2= @r1iB;n (t)i;n (t)@xB;n (t)xB;n (t) 0 x1 (t)= p(xB;n (t) 0 x1 (t))2 + (yB;n (t) 0 y1 (t))2 + (zB;n (t) 0 z1 (t))2
xB;n (t) 0 xi (t)0p2(xB;n (t) 0 xi (t)) + (yB;n (t) 0 yi (t))2 + (zB;n (t) 0 zi (t))2= @r1iB;n (t)i;n (t)@yB;n (t)yB;n (t) 0 y1 (t)= p(xB;n (t) 0 x1 (t))2 + (yB;n (t) 0 y1 (t))2 + (zB;n (t) 0 z1 (t))2yB;n (t) 0 yi (t)0p2(xB;n (t) 0 xi (t)) + (yB;n (t) 0 yi (t))2 + (zB;n (t) 0 zi (t))2= @r1iB;n (t)i;n (t)@zB;n (t)zB;n (t) 0 z1 (t)= p(xB;n (t) 0 x1 (t))2 + (yB;n (t) 0 y1 (t))2 + (zB;n (t) 0 z1 (t))2zB;n (t) 0 zi (t)0p2(xB;n (t) 0 xi (t)) + (yB;n (t) 0 yi (t))2 + (zB;n (t) 0 zi (t))2i;n(1.14)

これらは、適当に与えた利用者局の初期位置、および最小自乗法の逐次計算の途中段階における利用者
局の位置に対して、x, y , z 成分別に、それぞれの変化が距離差にどの程度の影響を及ぼすかを示すも
のである。

8

そして、行列

A (t) を次のように定める。
0 (0)
2 (0)
BB 2
BB 3 (0) 3 (0)BB
BB 4 (0) 4 (0)
BB .
..
BB ..
.
BB (0)(0)
A (t) = B
BB (t)2 (t)BB 2
BB (t)3 (t)BB 3BB (t)
4 (t)
BB 4
..
BB ...
.
@
n

c
f

1

;n

;n

2;n

(0)

;n

;n

3;n

(0) 0

;n

;n

4;n

s;n

s;n

;n

;n

(0) 0 0
..
..
.
.
0 111
s;n (0)c
02;n (t)

;n

;n

3;n t

;n

;n

4;n t

n

s;n

(t)

s;n

f

() 0

0 111 0
C
.. C
0
. C
CC
c . . . .. C
. C
CC
f
C
... ...
0 C
CCC
1 1 1 0 fc CC
0 111 0 C
C
.. C
C
0
. C
CC
c . . . .. C
. C
CC
f
C
... ...
0 C
A
111 0 c

0
c
f

c
f

() 0 0
..
..
.
.
0 111
s;n (t)(t)(1.15)f

さらに、最小自乗法により次のような行列計算を行うと、それぞれのパラメータの修正量が計算できる。

0 1x (t) 1BBCBB 1y (t) CCC
BB 1z (t) CC
BB
CC
B
0A (t) A (t)101 A (t) B (t) = BB 1N12 (t) CCC
BB 1N13 (t) CCBB
CBB 1N14 (t) CCC
BB
CC
..
B@
CA
.
B;n

B;n

B;n

n

T

n

n

AB;n

T

n

AB;n

(1.16)

AB;n1N1sAB;n (t)式 (1.6) から式 (1.16) までを、例えば 1xB;n  0:01 cm, 1yB;n  0:01 cm, 1zB;n  0:01 cm となるまで繰り返すことにより、1 エポックの測位計算が完了する。以上の逐次計算で求められるものは、
xyz 三次元直交座標系における利用者局位置 (xB , yB , zB ) である。したがって、さらにこれを日常私達が使用している、地球上の緯度, 経度、および高さの座標系に変換する必要がある。しかし、ここではこの問題には触れない。

1.4 アンビギュイティ
次に、アンビギュイティについて詳述する。
真のアンビギュイティは搬送波波長の整数倍となり、整数値バイアスとも呼ばれる。しかし、前述の計算で求めた値は整数とはならず、エポック毎に値の変化する実数となるので、この解を FLOAT 解という。FLOAT 解は整数値アンビギュイティを求める過程において、アンビギュイティを実数値で測位計算を行うもので、測位結果にはアンビギュイティの決定誤差が含まれる。

図 1.6 に前述の測位計算アルゴリズムを使用して計算を行い、同時に求めたアンビギュイティの一例

を時系列で示す。なお、二重位相差の計算に使用した衛星間の直線距離が最も近い場合を (a) 、最も遠い場合 (b) としてそれぞれ示す。それぞれ時間の経過と共に、ある整数値に収束してゆく様子が分かる。

9収束するまでに要する時間は異なったが、他の二重位相差のアンビギュイティも、同様にある整数値に92

アンビギュイティ（波長）アンビギュイティ（波長）収束した。

衛星間距離: 約12,233 km

90
88
86
84
0

10

20
30
40
時 間（分）

50

60

394
392
390
388
386
0

(a) 近距離衛星使用

衛星間距離: 約42,435 km

10

20
時

30
40
間（分）

50

60

(b) 遠距離衛星使用

図 1.6: アンビギュイティの収束状況
図 1.7 に二重位相差の計算に使用した 2 衛星間の直線距離と、アンビギュイティの標準偏差との関係を示す。同図から、衛星間の直線距離が短いほど、標準偏差が小さいことが分かる。この原因は、使用している衛星同士の距離が近いほど、電波の伝送路中の電離層、および対流圏の影響が同程度であり、アンビギュイティσ（波長）

二重位相差によって相殺される割合が高いためであると考えられる。

1.2
1.0
0.8
0.6
0.4
0.2
0.0
0

10000 20000 30000 40000
衛星間距離（km）

図 1.7: アンビギュイティの標準偏差と衛星間距離

FLOAT 解のままでは、測位結果にアンビギュイティの決定誤差が含まれる。そこで、求めたこの値を整数化する必要がある。整数化にあたっては、最も簡単な方法は四捨五入することである。しかし、実数値のアンビギュイティの標準偏差が大きい場合には、この方法は適当ではない。一般的には、アンビギュイティの平均値と、標準偏差から判断して、可能性のある整数値の、幾つかの組み合わせを仮定する。次に、最小自乗法の残差を計算して、その残差を最も小さくする組み合わせを、真の整数値のア
ンビギュイティとする方法が考えられる。
サイクル スリップが発生しなければ、アンビギュイティの値は変化しないという性質がある。これ

10

から、真のアンビギュイティが決定できれば、以後アンビギュイティは未知のパラメータとはせずに固定できることが分かる。したがって、エポック毎に式 (1.6) のように独立な方程式（独立な二重位相差）

が定義できる。一方、未知数は全部で 3 つ（未知の利用者局の xyz 三次元位置座標）なので、単一のエポックのみのデータを使用するだけで、解を得ることができる。さらに、アンビギュイティを整数値に固定することによって、アンビギュイティの決定誤差が排除できる。したがって、前述の搬送波位相デー
タによる測距精度に DOP を乗じた測位精度が期待できることになる。こうして得られた RTK-GPS の測位解を FIX 解といい、この呼び方は広く使用されている。

1.5 RTK-GPS 測位計算アルゴリズム - FIX 解
次に、RTK-GPS の FIX 解についての測位計算アルゴリズムを記述する。
各衛星からの搬送波の周波数が等しい（ド リフト無し ）と仮定する。すると前出の二重位相差の式(1.4) は、エポック t において式 (1.17) のように書き表すことができることはすでに記述した。
ccDDijAB (t) = (rjB (t) 0 rjA (t)) 0 (riB (t) 0 riA (t)) + NijABff= (rjB (t) 0 riB (t)) 0 (rjA (t) 0 riA (t)) + c NijAB
f(1.17)

ここで、前述の FLOAT 解の場合との差異は、アンビギュイティが既知で、しかも時間に無関係な、ある固定整数値であることである。

FLOAT 解において、計算した各二重位相差毎の整数値アンビギュイティは既知であり、基準局の位置も既知であると仮定して、未知のパラメータを整理する。衛星位置 (xi (t), yi (t), zi (t)), (xj (t), yj (t),zj (t)) は衛星から送信される航法メッセージによって既知である。したがって、未知のパラメータは利
用者局の xyz 三次元直交座標 (xB (t), yB (t), zB (t)) ということになる。式 (1.5) の右辺の整数値アンビ

ギュイティは、衛星からの電波を受信し続けていれば変化せず、二重位相差毎に一定の整数値である。
したがって、独立な二重位相差から定義される方程式の数よりも、未知のパラメータの方が少ないので、単一のエポックで利用者局の位置を求めることができる。

まず、式 (1.4) 、および (1.5) を使用して、衛星 1 を基準衛星とし、測位に使用した衛星数 s より 1 つ

少ない数の独立な二重位相差、すなわちエポック t (定義される。

8 c>> (DD12>> f>> c (DD13>< f
c>> f (DD14>>>>>: c (DD1f 1) において、独立な方程式が式 (1.18) のように

AB(t) 0 N12AB ) + (r2A (t) 0 r1A (t)) = r2B (t) 0 r1B (t)AB(t) 0 N13AB ) + (r3A (t) 0 r1A (t)) = r3B (t) 0 r1B (t)AB(t) 0 N14AB ) + (r4A (t) 0 r1A (t)) = r4B (t) 0 r1B (t)

sAB

..
.

..
.

(1.18)..
.

(t) 0 N1sAB ) + (rsA (t) 0 r1A (t)) = rsB (t) 0 r1B (t)しかし、FLOAT 解のときと同様に、線形方程式にようには解くことはできない。そこで、まずエポッ
ク t における n 回目の逐次演算後の利用者局位置 (xB;n (t), yB;n (t), zB;n (t)) とする。そのときの衛星 1（基準衛星）と利用者局 B 間、および衛星 s と利用者局 B 間の近似距離差 r1sB;n (t) = rsB;n (t)0 r1 (t)B;nは、式 (1.9) で表される。利用者局位置は初めに適当に与えるので、式 (1.6) 、および (1.18) は実際には両辺に対して等号は成立しない。そこで、これらの各式において（左辺 ）-（右辺）を計算して、次の11

ような行列B (t) を次のように定める。
0
BB fc (DD12 (t) 0 N12
BB 
BB c (DD13 (t) 0 N13
BB  f
B (t) = B
BB c (DD14 (t) 0 N14
BB f
..
BB
.
B@ 
c
n

AB

AB;n

AB

AB;n

AB

AB;n

n

f

1
) + r12 (t) 0 r12 (t) CCCC) + r13 (t) 0 r13 (t) CCCCC) + r14 (t) 0 r14 (t) C
CC
CC
..
..
.
.
CC

A


A

B;n

A

B;n

A

B;n

(DD1sAB (t) 0 N1sAB;n ) + r1sA (t)

0 r1

sB;n

(1.19)(t)

また、利用者局の近似位置 (xB;n (t), yB;n (t), zB;n (t)) に対する修正値を、それぞれ 1xB;n (t), 1yB;n (t),1zB;n (t) とすると、両者の関係は式 (1.12) のようになる。
次に、式 (1.9) から概略の利用者局位置から見た、衛星へ向かう視線の方向余弦の差の、各 xyz 座標@r1sB;n (t) @r1sB;n (t) @r1sB;n (t),,を計算すると、次式のようになる。
成分である@xB;n (t) @yB;n (t) @zB;n (t)

8
>>
>>
>>
>>
>>
>>
>>
>>
><
>>
>>
>>
>>
>>
>>
>>
>>
>:

(t) = @r1sB;n (t)@xB;n (t)xB;n (t) 0 x1 (t)= p
2(xB;n (t) 0 x1 (t)) + (yB;n (t) 0 y1 (t))2 + (zB;n (t) 0 z1 (t))2
xB;n (t) 0 xs (t)0p(xB;n (t) 0 xs (t))2 + (yB;n (t) 0 ys (t))2 + (zB;n (t) 0 zs (t))2= @r1sB;n (t)s;n (t)@yB;n (t)yB;n (t) 0 y1 (t)= p
2(xB;n (t) 0 x1 (t)) + (yB;n (t) 0 y1 (t))2 + (zB;n (t) 0 z1 (t))2yB;n (t) 0 ys (t)0p(xB;n (t) 0 xs (t))2 + (yB;n (t) 0 ys (t))2 + (zB;n (t) 0 zs (t))2= @r1sB;n (t)s;n (t)@zB;n (t)zB;n (t) 0 z1 (t)= p
2(xB;n (t) 0 x1 (t)) + (yB;n (t) 0 y1 (t))2 + (zB;n (t) 0 z1 (t))2zB;n (t) 0 zs (t)0p(xB;n (t) 0 xs (t))2 + (yB;n (t) 0 ys (t))2 + (zB;n (t) 0 zs (t))2s;nそして、行列A (t) を次のように定める。
0 (t)
BB 2
BB 3 (t)BB 4 (t)A (t) = B
BB .
B@ ..
n

;n

;n

n

;n

s;n

(t)

()
3;n (t)4;n (t)..
.
s;n (t)2;n t()3;n (t)4;n (t)..
.
s;n (t)2;n t(1.20)

1

1:0
C
1:0 C
CC
CC
1:0 C
.. C
C
. C
A
1:0

(1.21)

さらに、最小自乗法により、次のような行列計算を行う。すると、それぞれのパラメータの修正量が

12

0 1x
B
0A (t) A (t)101 A (t) B (t) = BBB 1yBB 1z@計算できる。

1
CC
CC
CC
A

(t)
B;n (t)B;n (t)
B;n

n

T

n

n

T

n

(1.22)

1T imeB;n (t)式 (1.22) までを、FLOAT 解の場合と同様に、例えば 1xB;n  0:01 cm, 1yB;n  0:01 cm, 1zB;n 
0:01 cm となるまで繰り返す。以上の過程により、1 エポックの測位計算が完了する。

1.6 KGPS/RTK-GPS 実例

10
0
–10
–20

10
0
–10
–20

–20

–10
0
10
経 度 方 向 (cm)

20

–20

(a-1) 水平方向分布

σ = 1.69 cm

7衛星
59.8
さ (m)

59.8
59.6
59.4

高

さ (m)

20

60
7衛星

高

–10
0
10
経 度 方 向 (cm)(b-1) 水平方向分布

60

59.2
59
0

2drms = 2.07 cm

7衛星
20 5秒毎

7衛星
5秒毎

緯 度 方 向 (cm)

緯 度 方 向 (cm)

20

59.6
59.4
59.2

10

20
時

30
40
間（分）

50

59
0

60

10

20
時

30
40
間（分）

(a-2) 高さの時系列(b-2) 高さの時系列(a) FLOAT 解(b) FIX 解

5060

図 1.8: KGPS 測位結果
ここで、実際に計算した FLOAT 解と、FIX 解の測位精度について考察する。図 1.8 (a), (b) に FLOAT

解、および FIX 解のうち水平方向分布と、高さの時系列の一例をそれぞれ示す。同図はまったく同じ搬

13送波位相データを使用して計算したものである。しかし、明かに FLOAT 解より FIX 解の精度が良い

ことが分かる。また、FLOAT 解は測位計算の開始時には特に精度が悪く、時間の経過と共にある一定の値に収束することが分かる。これは、衛星の移動によって利用者局の位置における「位置の面」の、見掛けの交角が変化し、搬送波位相データの測距誤差の、測位計算への影響が小さくなったためであると考えられる。

次に、RTK-GPS の測位精度の性能限界を調査するために、ゼロ ベース ライン6 （以下 ZBL と記述する）と基線長 20 m で測位を行った。図 1.9(a), (b) にこのときの測位結果の水平方向分布と、高さ方向の時系列の一例をそれぞれ示す。使用した GPS 受信機は、Trimble Navigation（以下トリンブルと記述する）社製の 4000SSi（以下 Trimble 4000SSi と記述する）である。この測位では、基準局から利用者局への測位用データの伝送は、受信機の入出力ポート同士をケーブルで直結した、最も理想的な状況で行っている。ZBL では、マルチパスが基準局と利用者局の搬送波位相データに全く同じ影響を及ぼすので、二重位相差の計算によって相殺され、マルチパスの影響の無い測位が可能となる。これに対（ ZBL の精度）+（マルチパス誤差）の精度で測位される。
して、基線長 20 m の場合では、

32

測位時間: 8.5時間
FIX率: 100.0％
2drms = 0.116 cm
ゼロ・ベースライン

059.28 ゼロ・ベースライン: σ = 0.084 cm

(m)

–1

–3
–4

さ

–2
基線長: 20 m
測位時間: 7.7時間
FIX率: 99.8％
2drms = 0.880 cm–5–5 –4 –3 –2 –1

0

1

59.26
59.24

高

緯 度 方 向 (cm)

1

59.22

2

基線長20 m: σ = 0.614 cm

3

59.20
0

10

経 度 方 向 (cm)(a) 水平方向分布

20

30
40
時 間 (分)5060(b) 高さの時系列

図 1.9: RTK-GPS 測位結果（ ZBL および基線長 20 m ）表 1.1 にこの測位結果をまとめて示す。ここで、FIX 率とは毎秒の正秒の全測位時間（ 3,600 回/時間）に対する、FIX 解の割合である。2drms と高さの標準偏差は共に mm オーダであり、これが Trimble

4000SSi の測位の性能限界であるといえる。

表 1.1: RTK-GPS 測位結果（ ZBL ）

6

基線長

測位時間

ZBL
20 m

8.5 時間7.7 時間

FIX 率

2drms

100.0 ％ 0.116 cm
99.8 ％ 0.880 cm

高さ標準偏差0.084 cm
0.614 cm

Zero base line（ ZBL ）： 基準局と利用者局で同一のアンテナを使用して、基線長を 0 としたもの。

14

1.7 RTK-GPS の幾何学的解釈
RTK-GPS の幾何学的にな意味合いについて考察する。

図 1.10 に示すように、RTK-GPS は二重位相差の計算に使用した 2 衛星を焦点とした回転双曲面の、複数組の交点として位置を求める方法である。そこで、まず衛星 1 (x1 , y1 , z1 ), 衛星 2 (x2 , y2 , z2 ) を焦点とする回転双曲面の式を求めることからアプローチを試みる。この 2 衛星と、ある点 (x, y , z ) との距離差を r とすると、この距離差が一定であるような点の集合を表す図形の式は、

F (x; y; z )pp= (x 0 x1 )2 + (y 0 y1 )2 + (z 0 z1 )2 0 (x 0 x2 )2 + (y 0 y2 )2 + (z 0 z2 )2 0 r= 0(1.23)となる。この式によって定義される図形は、2 衛星（回転双曲面の焦点となっている）を通る平面で切断した切り口にできる双曲線を、2 衛星を通る直線を回転軸として、その周りに回転させてできる回転面となる。これが回転双曲面体（ Hyperboloid ）と呼ばれる理由である (3) 。

図 1.10: 回転双曲面群の移動

1.8 RTK-GPS のデータ伝送フォーマット
RTK-GPS のデータ伝送手段としては、特定小電力無線システムがよく使用される。しかし、通信費用が高価ながら携帯電話回線も使用される。また、最近になってデ ィジタル MCA（ Multi-channel
Access ）業務用移動体無線通信システムを使用したものも実用化されている。しかし、依然としてデータ伝送システムのメディアとして、伝送に使用可能な周波数割り当てが整っているとはいえない状況にある。したがって、データ伝送フォーマットを作成する際には、データ量を小さくすることが望ましい。
このデータ伝送フォーマットに関しては、現在まで受信機メーカ各社共に独自のものを採用しており、公開されていない状況であった。このように、公表された標準フォーマットが無い状態では、利用者は15

異なるメーカの受信機を混在して使用することができないばかりか、測位結果に対して信頼性を確保することが困難な場合も多い。

RTK-GPS の測位用データのために、、RTCM SC-1047 がバージョン 2.1 文書中のメッセージ タイ
プ 18〜21 で、伝送フォーマットを定義している (4) 。受信機が RTCM の規定する RTK-GPS データに対応していれば、異なるメーカの受信機を使用しても RTK-GPS が可能である。一方、トリンブルは
ION GPS-96 で、伝送フォーマット CMR8 を公開した (5) 。以上のように、異なるメーカの受信機を混在して使用しても、RTK-GPS 測位ができる環境が整備されつつある。データ伝送フォーマットを規格化することの利点は、利用者としても測位用のデータの内容を把握した上で使用することにより、信頼性が向上することである。

(1) RTCM SC-104 Ver.2.1
DGPS を手軽るに利用するために、データ内容および伝送フォーマット等の規格が SC-104 Ver.2.1
で細かく決められており、現在このフォーマットはほぼ世界標準となっている。RTK-GPS に関するデータも Type18〜21 に含まれており、すでにこのフォーマットに準拠した GPS 受信機が市販されている。このデータを使用した RTK-GPS を RTCM RTK-GPS と呼ぶことにする。RTCM
RTK-GPS では Type18, 19 に搬送波位相データと、擬似距離測定データの生データの伝送フォー
マット等が定義されている。また、データ伝送が不安定な場合、RTK-GPS の測位結果がが取得できない可能性があることから、通常の DGPS 補正データ Type1 、および Type3 も同時に伝送される場合がある。この結果、RTCM RTK-GPS 測位では、データ伝送の状況が悪い場合にも測位は中断せず、自動的に DGPS 測位に切り替わる設定となっている場合が多い。
表 1.2 に二周波で Type18 を使用した場合の、衛星数とデータ量の関係を示す。Type19〜21 もまったく同様のデータ量である。ここでの 1 ワードは 30 ビットである。今後、公表される Ver. 3.0の文書では、実用化に向けたデータ伝送量の少ないフォーマットが公開されると聞いており、世界標準となる可能性が高いと考える。

表 1.2: 衛星数と RTCM Type18 のデータ量
衛星数

ワード 数 (words)

ビット 数 (bits)

4機

2+1+4 × 4 機=

19

19 × 30

=

570

5機

2+1+4 × 5 機=

23

23 × 30

=

690

6機

2+1+4 × 6 機=

27

27 × 30

=

810

7機

2+1+4 × 7 機=

31

31 × 30

=

930

8機

2+1+4 × 8 機=

35

35 × 30

=

1,050

9機

2+1+4 × 9 機=

39

39 × 30

=

1,170

10 機

2+1+4 × 10 機=

43

43 × 30

=

1,290

11 機

2+1+4 × 11 機=

47

47 × 30

=

1,410

12 機

2+1+4 × 12 機=

51

51 × 30=

1,530(2) Trimble Navigation CMR
前述の通り、RTCM では、RTK-GPS 測位用データの伝送フォーマットを Type18〜21 に定義しているが、データ量が多いため、9,600 bps 以上の伝送速度を実現できるデータ伝送システムを準備する必要がある。したがって、現在までのところ、RTK-GPS のデータ伝送フォーマットは、各メーカが独自に開発したものを使用しており、初期の段階ではどのメーカも、そのフォーマット

を公表していなかった。しかし、1996 年 9 月、カンザスシティで開催された ION GPS-96 において、ト リンブル社（ Trimble Navigation Co.,Ltd. ）は独自の RTK-GPS 測位用データ伝送フォー

マットである CMR を公開した。このフォーマットには、メッセージ プロトコルだけでなく、デー

7Radio Technical Commission for Maritime Services Special Committee No.104：米国海上無線技術委員会 第 104 特別会議
8 Compact measurement record： コンパクト観測レコード

16

タの圧縮伸長アルゴリズムが含まれている。これにより、RTCM RTK-GPS データと比較して半分以下のデータ伝送速度（ 2,400 bps ）で RTK-GPS 測位を実現できるといわれている。

CMR フォーマットの Type1 は、基準局のアンテナ位相中心の ECEF 座標と、そのオフセット量を含んでいる。基準局の状態を記述するパラメータには、基準局も移動することを前提にした「 Kinematic 」のパラメータが含まれている。つまり、そのような状態での RTK-GPS 測位が可能であることを示唆している。また、GPS 受信機の種類のパラメータに、他社の受信機の種類を記述するスペースがもうけられており、異なるメーカの受信機による RTK-GPS に対応する準備があることを示している (5) 。

同一のアンテナに接続した 4000SSi から、同時にそれぞれ 3 時間 25 分収集した平均データ量は 、RTCM RTK-GPS データ（ Type3, 18, 19, 59 ）では毎秒およそ 3.9 kbit/秒（ 1 セット ）であった。これに対して、CMR では約 1/4 の 1.0 kbit/秒（ 1 セット ）程度であった。したがって、CMR に要求される回線のデータ伝送性能は、RTCM RTK-GPS のそれと比較して、低速度で良いことになる。

17

参考文献
(1) 日本測地学会： 「新訂版 - 人工衛星による精密測位システム - 」, 社団法人 日本測量協会, 1989 年
11 月 15 日(2) 土屋 淳・辻 宏道： 「 GPS 測量の基礎」, 社団法人 日本測量協会(3) 占部 実・光藤 富士男： 「理工系一般教育 代数・幾何教科書」, 共立出版株式会社, 昭和 44 年 3 月
15 日(4) Radio Technical Commission for Maritime Services Special Committee No.104, ' RTCM Recommended Standards for Di erential Navstar GPS Service Version 2.1 ', RTCM Paper194-93/SC104STD, Washington,D.C., JANUARY 3 (1994)(5) Dr.Nicholas, C.Talbot, 'Compact Data Transmission Standard for High-Precision GPS', PROCEEDING OF THE 9TH INTERNATIONAL TECHNICAL MEETING OF THE SATELLITEDIVISION OF THE INSTITUTE OF NAVIGATION ION GPS-96 PART 1 OF 2, pp.861-871