リアルタイムキネマティック(RTK)測位の基礎と実装 Ba

リアルタイムキネマティック
(RTK)測位の基礎と実装
Basics and Implementations of Real-Time
Kinematic (RTK) Positioning
技術コンサルタント高須知二
Tomoji TAKASU
内容
• RTK測位の基礎
• 整数不定性決定
• RTK測位アルゴリズムの実装
http://gpspp.sakura.ne.jp
RTK測位
• リアルタイム移動体精密測位
• 搬送波位相観測値を使った相対測位
→測位解：相対座標(基線ベクトル)
• 移動観測点(Rover)+
基準観測点(Reference Station)
• データ通信リンク
• 高速整数不定性決定(OTF)
http://gpspp.sakura.ne.jp
RTK測位の概念
GPS/GNSS
Satellites
Rover
Baseline
Reference
Station
Communication Link
http://gpspp.sakura.ne.jp
RTK測位の基礎
• 擬似距離観測値
• 搬送波位相観測値
• 航法メッセージ・放送暦
• 幾何学距離
• 単独測位
• 二重位相差
• 最小二乗法・カルマンフィルタ
http://gpspp.sakura.ne.jp
擬似距離測定値
測位コードにより測定された信号伝搬時間
×光速
Pi  c(t r  t )   Pi 
s
  c(dt  dT )  I i  T   Pi
測位コード(PRNコード)
tr
http://gpspp.sakura.ne.jp
搬送波位相測定値
GPS/GNSS搬送波位相角の連続測定値
Li  iΦi 
  c(dt  dT )  I i  T  i N i   Li
搬送波位相バイアス
s
Ni  0r ,i  0,i
 ni （整数不定性）
http://gpspp.sakura.ne.jp
航法メッセージ・放送暦
• 測位信号に重ね合わされて送信
• 概ね2時間間隔で更新
• 概略衛星軌道(Almanac)
• 衛星軌道(Ephemeris:放送暦)
• 衛星時計誤差
• 衛星ヘルス情報,電離層情報 etc
• 基準時刻：GPS Week #+Toe, Toc
http://gpspp.sakura.ne.jp
衛星軌道(放送暦)(1/3)
Satellite
tk
0:00
GPS Week
day 1
t oe1
t oe 2 t
t oe 3 GPST(sec)
http://gpspp.sakura.ne.jp
衛星軌道(放送暦)(2/3)
Satellite
E
A
Ae
Satellite
z
 e
rs
r
u
u

y

x
i
http://gpspp.sakura.ne.jp
衛星軌道(放送暦)(3/3)
t k  t  t oe, n   / A3  n, M  M 0  ntk
E  M  e sin E (ケプラー方程式)
  ATAN2( 1  e 2 sin E , cos E  e)  

u  
  Cus Cuc 
  
 
 sin 2 

 r    A(1  e cos E )    C rs C rc 
cos 2 

 i   i  IDOT t   C

Cic 
k   is
   0
   0  (   e )t k   e t oe
r s (t )  Rz ( ) Rx (i )(r cos u, r sinu, 0)T
http://gpspp.sakura.ne.jp
衛星時計誤差
dT (t )  a f 0  a f 1 (t  toc )  a f 2 (t  toc )
s
2
 t r  tGD
相対論補正
t r 
 2  Ae sin E
擬似距離バイアス補正
  TGD

t GD    TGD
 0

c2
( L1)
2
2
( L 2) (  f1 / f 2 )
( LC )
http://gpspp.sakura.ne.jp
幾何学距離
光差方程式(Light-Time Eq.)
s
r

Satellite s
s
rr
rrs  rr  RZ ( e  rs / c) r s (t   rs / c)
(t  t r  dtr )
( rs / rrs )T  e rs  rrs /  rs
Elrs   arcsin(e rs  rr / rr )
e rs
 rs
s
Elr
Receiver r
http://gpspp.sakura.ne.jp
最小二乗法
観測方程式(非線形)
y  h( x)  ε  h( x0 )  H ( x  x0 )  ε
( H  h( x) / x x  x , E{ε}  0, E{εεT }  Q y )
0
最小二乗推定値
ˆx  x0  ( H T Q y 1H ) 1 H T Q y 1 ( y  h( x0 ))
推定値共分散行列
Qxˆ  ( H
T
1
1
Qy H )
http://gpspp.sakura.ne.jp
カルマンフィルタ
逐次最小二乗法
推定値
観測データ
y1 , y 2 ,...
観測更新
xˆ1 , xˆ 2 ,...
時間更新
http://gpspp.sakura.ne.jp
単独測位(1/3)
Satellite 2
Satellite 3
Satellite n
Satellite 1
Pr2
Pr1
Pr3
Prn
Receiver r
http://gpspp.sakura.ne.jp
単独測位(2/3)
未知パラメータ
x  (rr , cdtr )
T
観測モデル
T
 e1T
  1r  cdtr  cdT1  Tr1 
 r


 e 2T
  r2  cdtr  cdT 2  Tr2 
r

 3

h( x )    r  cdtr  cdT 3  Tr3  H   e 3T
r




 
 n

  r  cdtr  cdT n  Trn 
 e nT


 r
1

1
1


1
http://gpspp.sakura.ne.jp
単独測位(3/3)
観測量
y
1
2
3
n T
 ( Pr , Pr , Pr ,..., Pr )
( Prs  C1Prs1  C2 Prs 2 )
最小二乗解(Gauss-Newton)
xˆ  (0,0,0,0)
0
T
ˆx i 1  xˆ i  ( H T H ) 1 H T ( y  h( xˆ i ))  xˆ
http://gpspp.sakura.ne.jp
単独測位解例
East (m)
10
Single Point Positioning Solution : 960583
MEAN: 0.8191m RMS: 1.4907m
0
North (m)
-10
10
MEAN: 1.0867m RMS: 1.4880m
0
Up (m)
-10
10
MEAN: 0.4081m RMS: 2.7906m
0
-10
0:00
1:00
2:00
3:00
http://gpspp.sakura.ne.jp
二重位相差
AB
Lab
A
 ( La
Satellite A
LaA
Receiver a
B
A
 La )  ( Lb
LbA
LBa
B
 Lb )
Satellite B
LBb
Receiver b
http://gpspp.sakura.ne.jp
二重位相差観測モデル
AB
Lab
AB
  ab
AB
  ab
AB
 c (dtab
AB
AB
 dTab )  I ab
AB
 Tab
AB
 N ab
AB
dtab
 (dtaA  dtaB )  (dtbA  dtbB )  0
dTabAB  (dTaA  dTaB )  (dTbA  dTbB )  0
AB
AB
N ab
 ( N aA  N aB )  ( NbA  NbB )  nab
AB
Lab
AB
AB
AB
AB
AB
  ab  I ab  Tab  nab   ab
AB
AB
AB
AB
Lab   ab  nab   ab （短基
線）
http://gpspp.sakura.ne.jp
RTK測位
Satellite 2
Satellite 1
Receiver r
(Rover)
Satellite s
(Reference Sat)
Satellite n
Receiver R
(Reference
Station)
http://gpspp.sakura.ne.jp
RTK測位(1/5)
未知パラメータ
T
T T
s1
s2
sn T
x  (rr , N ) , N  ( N rR , N rR ,..., N rR )
(観測点座標+搬送波位相バイアス二重差)
未知パラメータ/共分散初期値
s1
s2
sn
xˆ1  (rPP (t1 )T , N rR
(t1 ), N rR
(t1 ),..., N rR
(t1 ))T
sk
sk
sk
(rPP (t1 ) : t1　単独測位解 N rR (t1 )  ( LrR (t1 )  PrR (t1 )) /  )

P1  diag( r ,  r ,  r ,  N ,  N ,...,  N )
2
2
2
2
2
2
http://gpspp.sakura.ne.jp
RTK測位(2/5)
二重位相差観測モデル(短基線)
T
s1
s1 
s

  rR


N
rR 
 er

s2
s2
 e sT
  rR

 N rR
h( x )  
 H  r




 sT
  sn  N sn 
rR 
 rR
 er
1T
 er
2T
 er
 0  0 
0   0 


nT
 er
0
  

0  
http://gpspp.sakura.ne.jp
RTK測位(3/5)
二重位相差観測量/観測誤差共分散
 LsrR1 
 4 2



 LsrR2 
 2 2
yk  
 Rk  
  
 
 2 2
 Lsn 

 rR 
2 2  2 2 
4 2  2 2 



 
2 2  4 2 
http://gpspp.sakura.ne.jp
RTK測位(4/5)
拡張カルマンフィルタ観測更新


K k  Pk H k ( H k Pk H k  Rk )
T
T

1

xˆ k  xˆ k  K k ( y k  h( xˆ k ))
Pk  ( I  K k H k ) Pk

xˆ k : t k 推定値
Pk : t k 推定値共分散
http://gpspp.sakura.ne.jp
RTK測位(5/5)
拡張カルマンフィルタ時間更新
ˆx k 1  Φxˆ k  (rPP (t k 1 )T ,0,0,..., 0)T

Pk 1  ΦPk Φ T  Q
Φ  diag(0,0,0,1,1,...,1)
Q  diag( r ,  r ,  r ,0,0,..., 0)
2
2
2
http://gpspp.sakura.ne.jp
RTK測位解例(FLOAT解)
RTK Position (FLOAT) : Baseline=960583-92110 (14.1km)
East (m)
0.2
MEAN:-0.0163m RMS: 0.0423m
0
North (m)
-0.2
0.2
MEAN:-0.0022m RMS: 0.0485m
0
Up (m)
-0.2
0.2
MEAN:-0.0073m RMS: 0.0517m
0
-0.2
0:00
1:00
2:00
3:00
http://gpspp.sakura.ne.jp
整数不定性決定(AR)
• 搬送波位相バイアス二重差→整数
• 整数条件を利用し推定値を整数に固定す
ることにより一般に測位精度が向上する。
• 短基線、キネマティックで特に有効
• RTK:
初期化時間短縮、サイクルスリップ対応
→高速 (On-The-Fly:OTF)決定
http://gpspp.sakura.ne.jp
高速(OTF)AR手法の例
手法
開発者
LSAST
FARA
Hatch
Free, Beutler
Modified
Cholesky
Decomposition
Euler,
Landau
LAMBDA
Null Space
Teunissen
Martin-Neira
FASF
Chen,
Lachapelle
OMEGA
Kim, Langley
探索
方法
独立
一括
処理
エポック
単一
複数
探索空間
限定手法
なし
条件付
一括
複数
なし
一括
独立
複数
単一
変換/条件付
変換
一括
複数
条件付
独立単一/複数変換/条件付
D.Kim et al., GPS Ambiguity Resolution and Validation Methodologies, Trend and Issues, International Symposium on GPS/GNSS, 2000
http://gpspp.sakura.ne.jp
高速(OTF)AR手法
• 多数衛星の利用
• 多周波、擬似距離の利用
• 1エポックでの決定(瞬時決定)も可能
• 主な手法の分類
ワイドレーン法
整数最小二乗(Integer Least Square)
その他
http://gpspp.sakura.ne.jp
ワイドレーン法(1/3)
• 衛星ペア毎に整数不定性を独立に解く
• 近代化GPS/Galileoで有効
→3周波測位信号(L1/L2/L5)の利用
• ワイドレーン、エクストラワイドレーン
• 幾何学フリー観測モデル
http://gpspp.sakura.ne.jp
ワイドレーン法(2/3)
0  L1 / 1 
 L1 / 1   1 0

 


 LWL / WL    1  1 0  L2 / 2 
L
  0 1  1 L /  
/

 EWL EWL  
 5 5 
L1 (1  19.0cm) ：L1搬送波位相
LWL (WL  86.2cm) ：ワイドレーン
LEWL (EWL  5.86m) ：エクストラワイドレーン
http://gpspp.sakura.ne.jp
ワイドレーン法(3/3)
 P  P1
N EWL  round(( LEWL   P ) /  EWL )
 EWL  LEWL   EWL N EWL
NWL  round(( LWL   EWL ) / WL )
WL  LWL  WL NWL
N1  round(( L1  WL ) / 1 )
N 2  N1  NWL , N 5  N 2  N EWL
http://gpspp.sakura.ne.jp
整数最小二乗(1/3)
• 整数不定性を整数とした制約で最小二乗
条件を満たす解を探索
• 誤差=正規分布→最大正解率
• ステップ
– (実数)最小二乗解(FLOAT解)
– 整数不定性につき残差(実数解-整数解)二乗
和を最小とする整数解を探索
– 整数解を固定し再度最小二乗解(FIX解)
http://gpspp.sakura.ne.jp
整数最小二乗(2/3)
整数最小二乗条件
y  Aa  Bb  ε,
整数基線解
不定性 (実数)
(1)
min
n
aR ,bR
m
min
n
aZ ,bR m
y  Aa  Bb
最小二乗
2
2
Qy
y  Aa  Bb
 Qaˆ
ˆ
 aˆ , b, Q  
Q
FLOAT解  bˆaˆ
2
Qy
Qaˆbˆ 

Qbˆ 

 a 最近整数格子点探索
(2) minn aˆ  a Q
aˆ
aZ


1
ˆ
(3) b  b  Qbˆaˆ Qaˆ (aˆ  a ) FIX解

2
Q
 (  ) T Q 1 (  )
http://gpspp.sakura.ne.jp
整数最小二乗(3/3)

最近整数格子点探索 a  arg min aˆ  a
n
aZ
 0.4 

aˆ  
aˆ  a
  0.1
 1.6 2.8 

Q  
 2.8 5.6 
Qaˆ
2
Qaˆ
aˆ
1

a
50
: FLOAT解
: FIX解
2
0
0
-1
-1
0
a1
1
-2
a2
2
http://gpspp.sakura.ne.jp
LAMBDA
• Least-square AMBiguity
Decorrelation Adjustment
(Teunissen, 1995)
• 整数最小二乗を利用したAR手法の一つ
• Z変換(無相関化)による前処理
→探索効率化
• 最近整数格子点探索アルゴリズム
http://gpspp.sakura.ne.jp
Z変換(無相関化)
• 最近整数格子点探索の効率化
• 整数格子点⇔整数格子点変換
(unimodular変換)
• LAMBDA Decorrelation(無相関化)
整数Gauss変換+行入れ替え
n
n
T
a

Z

z

Z
zZ a
Q zˆ  Z Qaˆ Z
T
aˆ  a
2
Qaˆ
 zˆ  z
2
Q zˆ
http://gpspp.sakura.ne.jp
Z変換の例
  1 0  T   1 0 
, Z  

Z  
2  1
 2  1


z
a
T
1
Z
T
0
2
1
aˆ
zˆ
-1
Z
-2
T
0
-1
-1
0
1
2
-2
-1
0
1
http://gpspp.sakura.ne.jp
最近整数格子点探索(1/3)
2
④
1
②
0
zˆ

z
①
③
-1
-2
-1
0
1
http://gpspp.sakura.ne.jp
最近整数格子点探索(2/3)
Q zˆ  LT DL
( z  zˆ)T Q zˆ 1 ( z  zˆ)  ( z  zˆ)T L1 D 1 LT ( z  zˆ)   2
z  z  LT ( z  zˆ)
n
z n  zˆ n, zi  zˆi 
l
ji ( z j
 zˆ j ) (i  n  1, n  2,...,1)
j i 1
( z n  z n ) 2 ( z n 1  z n 1 ) 2
( z1  z1 ) 2

 ... 
 2
dn
d n 1
d1
http://gpspp.sakura.ne.jp
最近整数格子点探索(3/3)
-1
0
-1
0
-1
1
z4
0
z3
0
z2
1
0 -2 -1 0 1 2
z1
http://gpspp.sakura.ne.jp
最近整数格子点探索(3/3)
-1
0
-1
0
-1
1
z4
0
z3
0
z2
1
0 -2 -1 0 1 2
z1
http://gpspp.sakura.ne.jp
最近整数格子点探索(3/3)
-1
0
-1
0
-1
1
z4
0
z3
0
z2
1
0 -2 -1 0 1 2
z1
http://gpspp.sakura.ne.jp
最近整数格子点探索(3/3)
-1
0
-1
0
-1
1
z4
0
z3
0
z2
1
0 -2 -1 0 1 2
z1
http://gpspp.sakura.ne.jp
最近整数格子点探索(3/3)
-1
0
-1
0
-1
1
z4
0
z3
0
z2
1
0 -2 -1 0 1 2
z1
http://gpspp.sakura.ne.jp
LAMBDA手順
最小二乗解(FLOAT解)
aˆ , bˆ, Q
Z変換(無相関化)
zˆ  Z T aˆ, Qzˆ  Z T Qaˆ Z
最近整数格子点探索
逆Z変換

z  arg min zˆ  z
zZ n
2
Qzˆ

T 
aZ z


1
ˆ
ˆ
整数最小二乗解(FIX解) b  b  Qbˆaˆ Qaˆ (a  a)
http://gpspp.sakura.ne.jp
LAMBDA実行時間例
Execution Time (ms)
15
: with decorrelation
: without decorrelation
(Pentium 4 3.2GHz
Intel C/C++8.0)
10
5
0
5
10
15
20
25
30
N : Number of Ambiguities
35
40
http://gpspp.sakura.ne.jp
整数不定性の検定
• FIX解の統計的妥当性の検査
• 最大正解率/最小不正解(ミスFIX)率
• Ratio Test
 2
aˆ  a 2 Q


aˆ
次善解 ) FIX解 a 2 :  2   (a : aˆ  a Q
aˆ
• その他： F-ratio Test、W-Test etc
http://gpspp.sakura.ne.jp
RTK測位解例(FIX解)
East (m)
0.2
RTK Position (FIX) : Baseline=960583-92110(14.1km)
MEAN: 0.0137m RMS: 0.0158m
0
North (m)
-0.2
0.2
MEAN: 0.0031m RMS: 0.0089m
0
Up (m)
-0.2
0.2
MEAN:-0.0120m RMS: 0.0257m
0
-0.2
0:00
1:00
2:00
3:00
http://gpspp.sakura.ne.jp
RTK測位アルゴリズムの実装
• 多周波、擬似距離の利用
• 基準衛星の選択・切り替え
• サイクルスリップ・アウトライアの検出
• 電離層遅延推定
• 精密暦・精密補正の利用
• 複数基線の利用
• その他
http://gpspp.sakura.ne.jp
多周波、擬似距離の利用
• 複数周波数、擬似距離観測値の利用
→Ambiguity FIX時間の短縮
s1
s1
sn
sn T
x  (rr T , N1,2T )T , N1,2  ( NrR
,
N
,..,
N
,
N
rR 2
rR 1
rR 2 )
1
s1
s1 
 LsrR1 1 
  rR


N
1 rR 1 



s1
s1
 LsrR1 2 
  rR

 2 N rR
2

s1
y   PrRs11 , h( x )  
 rR

s1
 PrRs1 



rR
2





  


http://gpspp.sakura.ne.jp
基準衛星選択・切り替え(1/2)
基準衛星
Sat1
Sat2
Sat3
Sat4
Sat5
欠測
t
http://gpspp.sakura.ne.jp
基準衛星選択・切り替え(2/2)
• (仰角最大&欠測でない)衛星→基準衛星
• 推定値の切り替え(基準衛星: r→s)
x n  ( N r1 ,..., N rr ,..., N rs ,..., N rn )T , Pn
0  0 

0  0 
A  I 



0  0 

1  0

1  0



1  0 
x n '  Ax n  ( N s1 ,..., N sr ,..., N ss ,..., N sn )T , Pn '  APn AT
http://gpspp.sakura.ne.jp
サイクルスリップの検出
• アウトライア検出/除外：
事前残差/事後残差検査



ˆ
s  yk  h( x k ), s  yk  h( xˆ k )


T

Qs  H Pk H  Rk , Qs  Rk

2  
2 
si  n qii , si  n qii (n  3～5)
• 連続的なアウトライア
サイクルスリップ→推定値再初期化
http://gpspp.sakura.ne.jp
電離層遅延推定(1/4)
• 電離層遅延二重差の推定
→中基線(>20km)の測位精度向上
• 電離層遅延時間変動
→Random-Walkモデル
• 推定値の収束に時間がかかる
→初期化時間
http://gpspp.sakura.ne.jp
電離層遅延推定(2/4)
x  (rr , I
T
T
T T
, N1,2 ) , I
s1 s 2
sn T
 (I rR , I rR ,..., I rR )
s1
s1
s1
 LsrR1 1 




I


N
rR
rR
1
rR
1




s1
s1
2
2 s1
s1
 LrR 2 
  rR  f1 / f 2 I rR  2 N rR 2 

s1
s1
y   PrRs11 , h( x )  
 rR  I rR

s1
2
2 s1
 PrRs1 




f
/
f
rR
1
2 I rR
2





  


Φ  diag(0,0,0,1,1,...,1,1,1,...,1)
Q  diag( r 2 ,  r 2 ,  r 2 , t I2 , t I2 ,..., t I2 ,0,0,..., 0)
http://gpspp.sakura.ne.jp
電離層遅延推定(3/4)
0.2
推定 (FIX率=98.1%)
0固定(FIX率=100%)
RMS: 0.0307m
RMS: 0.0060m
RMS: 0.0183m
RMS: 0.0077m
RMS: 0.0370m
RMS: 0.0183m
0
-0.2
0.2
0
-0.2
0.2
0
-0.2
0:00
1:00
2:00
0:00
(Baseline=23.1km)
1:00
2:00
3:00
http://gpspp.sakura.ne.jp
DD of Ionospheric Delay (m)
電離層遅延推定(4/4)
0.1
0.05
0
-0.05
-0.1
0:00
1:00
(Baseline=23.1km)
2:00
3:00
http://gpspp.sakura.ne.jp
精密暦・精密補正の利用
• 精密暦：長基線(>100km)
IGS Ultra-Rapid(IGU)(予報値)
精度：<10cm、15分間隔値の補間
• 精密補正
対流圏遅延：高度差,長基線(>数10km)
アンテナPCV：異機種アンテナ間基線
局位置変動: 長基線(>数100km)
Phase-Windup：移動体,姿勢変動
http://gpspp.sakura.ne.jp
複数基線の利用
複数基線
単一基線
Rover
Reference
Station
http://gpspp.sakura.ne.jp
ネットワークRTK
• 基準局網による擬似基準局データ生成
• 基準局網による補正情報
電離層、対流圏、衛星軌道 etc.
• 主な方式
VRS (Trimble)
FKP (Geo++)
Multiref (Kalgary大)
http://gpspp.sakura.ne.jp
RTK用通信標準
• RTCM SC-104
DGPS補正情報
RTK用観測データ
ネットワークRTK用補正情報
etc.
• NTRIP(Networked Transport of
RTCM via Internet Protocol)
• NMEA
http://gpspp.sakura.ne.jp
まとめ
• RTK測位の基礎
擬似距離/搬送波位相、放送暦、最小二乗
幾何学距離、単独測位、RTK測位 etc
• 整数不定性決定
整数最小二乗、LAMBDA etc
• RTK測位アルゴリズムの実装
基準衛星選定、サイクルスリップ検出、
電離層推定 etc
http://gpspp.sakura.ne.jp
おまけ
n 2  1  f N2 / f 2 ( Appleton Hartee formula, L  band)
f N2 

N ee2
4  0 me
2
(electronplasmafreq.)
N e (1.602 1019 ) 2
12
4 8.854 10
2
31
 9.109 10
 80.60N e
n  1  f N2 / f 2  1  f N2 / 2 f 2    1  N I
N I  f N2 / 2 f 2    40.30N e / f 2
c refractavility)
(ionospheri
http://gpspp.sakura.ne.jp

Download Report