搬送波位相測定値による 精密測位の理論及び解析処理 Precise positioning

GPS/GNSSシンポジウム2005
搬送波位相測定値による
精密測位の理論及び解析処理
Precise positioning theory and analysis with
carrier-phase measurements
技術コンサルタント 高須 知二
Tomoji TAKASU
内容 Contents
• A.4.1 はじめに
• A.4.2 精密測位のモデル
• A.4.3 パラメータ推定の手法
• A.4.4 相対測位
• A.4.5 精密単独測位(PPP)
• A.4.6 精密解析の実際
• A.4.付 座標系の定義と変換
http://gpspp.sakura.ne.jp
A.4.1 はじめに
精密測位
Precise Positioning
GPS/GNSS搬送波位相測定値に
よる高精度測位技術
http://gpspp.sakura.ne.jp
A.4.2.1 擬似距離
測位コードにより測定された信号伝搬時間
×光速
P  c(t r  t )   P
s
観測時刻
(受信機時計)
送信時刻
(衛星時計)
測位コード(PRNコード)
tr
http://gpspp.sakura.ne.jp
A.4.2.1 擬似距離
測位コードにより測定された信号伝搬時間
×光速
P  c(t r  t )   P 
s
  c(dt  dT )  I  T   P
幾何学距離
電離層
観測誤差
衛星 遅延 対流圏 (10cm
受信機
遅延 ~数m)
時計誤差 時計誤差
http://gpspp.sakura.ne.jp
A.4.2.1 搬送波位相
GPS/GNSS搬送波位相角の連続測定値
L  Φ
搬送波波長(m) 搬送波位相(cycle)
搬送波
tr
http://gpspp.sakura.ne.jp
A.4.2.1 搬送波位相
GPS/GNSS搬送波位相角の連続測定値
L  Φ 
  c(dt  dT )  I  T  N   L
幾何学距離
電離層
衛星 遅延 対流圏
受信機
遅延
時計誤差 時計誤差
観測誤差
(1mm
~数cm)
http://gpspp.sakura.ne.jp
A.4.2.1 搬送波位相
GPS/GNSS搬送波位相角の連続測定値
L  Φ 
  c(dt  dT )  I  T  N   L
搬送波位相バイアス
N
s
 0r  0
 n (整数不定性)
http://gpspp.sakura.ne.jp
A.4.2.2 幾何学距離
y
Satellite
Receiver
z
x
r (t   )
s
送信時
http://gpspp.sakura.ne.jp
A.4.2.2 幾何学距離
y
Satellite
Receiver
rr (t )

z
x
RZ (E )r (t   )
s
 E
受信時
http://gpspp.sakura.ne.jp
A.4.2.2 幾何学距離
光差方程式(Light-Time Equation)を解く
  rr (t )  RZ ( E )r (t   )
s
観測点位置(受信時)
衛星位置(送信時)
   / c, t  t r  dt
電波伝搬時間
http://gpspp.sakura.ne.jp
A.4.2.2 幾何学距離
光差方程式(Light-Time Equation)を解く
  rr (t )  RZ ( E )r (t   )
s
   / c, t  t r  dt
受信時刻
受信機
観測時刻
(受信機時計) 時計誤差
http://gpspp.sakura.ne.jp
A.4.2.2 衛星方位・仰角
• 衛星方位・仰角 Az, El
z (Up) Satellite
r '  E r ( Rz ( E )r (t   )  rr (t ))
s
rs'
s
 ( x s , y s , z s )T
Az  ATAN2( x s , y s )
El  arcsin
z
s
Receiver
El
Az
rs'
x (East)
y (North)
http://gpspp.sakura.ne.jp
A.4.2.3 衛星軌道・衛星時計
r : 衛星位置
s
dT : 衛星時計誤差
• 放送暦 (Broadcast Ephemeris)
航法メッセージ中に含まれる
精度 : 軌道~数m, 時計~10nsec
• 精密暦 (Precise Ephemeris)
世界中の観測局データにより決定した
高精度衛星軌道・時計情報
http://gpspp.sakura.ne.jp
A.4.2.3 IGS精密暦
種類
Final
(最終暦)
決定/予報値
決定値
決定値
時間遅れ
13 日
17 時間 3 時間 即時
提供頻度
1 週間毎
1 日毎
6 時間毎
15 分
5分
<5cm
<0.1ns
15 分
5分
<5cm
0.1ns
15 分
15 分
<5cm 10cm
0.2ns 5ns
時間
間隔
精度
軌道
時計
軌道
時計
Rapid
Ultra-Rapid
(速報暦) (超速報暦)
決定
予報
http://gpspp.sakura.ne.jp
A.4.2.3 精密暦
• 衛星軌道
一定間隔の衛星位置(地球固定座標)
→多次多項式補間(10次以上)
• 衛星時計
一定間隔の衛星時計誤差(バイアス)
→補間誤差、時計の定義に注意
http://gpspp.sakura.ne.jp
A.4.2.4 電離層遅延
I
40.3
f2

Ndl 
40.3 10
16
f2
TEC
Total Electron Content
総電子数 (TECU=1016el/m2)
f : 搬送波周波数 (Hz)
http://gpspp.sakura.ne.jp
A.4.2.4 電離層フリーLC
二周波搬送波位相による電離層遅延項消去
LC  C1 L1  C 2 L2 
  c(dt  dT )  T  N LC   LC
C1 
f1
2
f1  f 2
2
2
, C2 
f2
2
f1  f 2
2
2
http://gpspp.sakura.ne.jp
A.4.2.4 電離層モデル
Single Layerモデル
Ionospheric
Pierce Point
z’
Satellite
1 40.3  1016
I
2
cos z '
f
 TEC(t ,  pp ,  pp )
H
z
Ionosphere
Receiver
Earth
RE
http://gpspp.sakura.ne.jp
GIM(Global Ionos Maps)
100
90
80
70
60
50
40
30
20
10
http://gpspp.sakura.ne.jp
0
A.4.2.5 対流圏遅延
静水圧遅延
湿潤遅延
T  M dry ( El)ZHD  M wet ( El)ZW D
天頂静水圧遅延
(2.0~2.4m)
天頂湿潤遅延
(0~0.4m)
M dry , M wet : マッピング関数
http://gpspp.sakura.ne.jp
A.4.2.5 対流圏遅延
• 天頂静水圧遅延(ZHD)
地上気圧(hPa)
ZHD 
0.0022768P0
1  0.00266cos 2  2.8  107 h
緯度
• 天頂湿潤遅延(ZWD)
未知パラメータとして推定
高度(m)
http://gpspp.sakura.ne.jp
A.4.2.5 マッピング関数
NMF(Niell Mapping Function)
1
M ( El) 
仰角 sin El 
a
Zenith
Satellite
b
1
1 c
a
b
sin El 
sin El  c
a,b,c : 係数(緯度、高度、通算日)
El
Receiver
http://gpspp.sakura.ne.jp
A.4.2.6 受信アンテナ位相中心
Satellite
Antenna
Phase Center
Geodetic
Reference
Point
Antenna
Phase Center
Offset
Antenna
Reference Point
(ARP)
http://gpspp.sakura.ne.jp
A.4.2.6 受信アンテナ位相中心
• アンテナ位相中心パラメータ
L1,L2 アンテナ位相中心オフセット
L1,L2 アンテナ位相中心変動(PCV)
• アンテナ機種毎テーブル
地上精密計測データベース
http://gpspp.sakura.ne.jp
A.4.2.7 衛星アンテナ位相中心
• 精密暦 : 衛星重心位置
• 衛星アンテナ位相中心オフセット補正
衛星
アンテナ位相中心(IGS)
x(m)
y(m)
z(m)
Block II/IIA
0.279
0.000
1.023
Block IIR
0.000
0.000
0.000
http://gpspp.sakura.ne.jp
A.4.2.7 衛星固定座標系
Sun
Satellite
y
x
z
Earth
http://gpspp.sakura.ne.jp
A.4.2.8 局位置変動
• 地球潮汐:
月・太陽の引力による弾性変形
固体地球潮汐(Solid Earth Tide)
海洋荷重(Ocean Loading)
極運動潮汐(Pole Tide) etc
• 局位置変動(Site Displacement)
最大上下変動:~20cm
http://gpspp.sakura.ne.jp
A.4.2.8 局位置変動
Up (m)
North (m)
East (m)
TSKB : 2004/1/1 0:00-1/2 0:00
0.2
0.1
0
-0.1
-0.2
0.2
0.1
0
-0.1
-0.2
0.2
0.1
0
-0.1
-0.2
1/1 0:00
1/1 6:00
1/1 12:00
1/1 18:00
1/2 0:00
http://gpspp.sakura.ne.jp
A.4.2.9 Phase Wind-up効果
• GPS搬送波
右旋円偏波
• 衛星ー受信機アンテナの相対的回転
→搬送波位相の進行、後退
http://gpspp.sakura.ne.jp
A.4.2.10 相対論効果
• 特殊相対論効果:衛星時計の遅れ
• 衛星軌道離心率 ≠ 0
→衛星時計誤差の周期変動
→最大で10m以上
• 放送暦・精密暦時計:相対論効果を除去
dT  dT0 
衛星時計誤差
2r  v
s
c
2
s
衛星位置・速度
(慣性座標系)
http://gpspp.sakura.ne.jp
A.4.2.11 マルチパス・観測誤差
• マルチパス
• 電離層遅延・対流圏遅延補正残差
• 受信機雑音
• カットオフ角
• 仰角によるWeighting
http://gpspp.sakura.ne.jp
A.4.2.12 線形結合
記号
線形結合
L1
L2
LC
LG
WL
NL
MW
MP1
MP2
(L1 搬送波位相)
(L2 搬送波位相)
電離層フリー
幾何学フリー
ワイドレーン
ナローレーン
Melbourne-Wübbena
L1 マルチパス
L2 マルチパス
波長
(cm)
19.0
24.4
86.2
10.7
86.2
-
電離
層
1.0
1.6
0.0
0.6
1.3
1.3
0.0
0.0
0.0
雑音
(cm)
0.3
0.3
0.9
0.4
1.7
1.7
21
30
30
http://gpspp.sakura.ne.jp
A.4.3.1 最小二乗法
観測方程式
z  h( x )  ε
観測量
観測誤差
観測モデル パラメータ
z  ( z1 , z2 , z3 ,..., zm )T
x  ( x1 , x2 , x3 ,..., xn )T
http://gpspp.sakura.ne.jp
A.4.3.1 最小二乗法
非線形最小二乗法
初期値 x0 の周りで線形化
z  h( x)  ε
 h( x0 )  H ( x  x0 )  ε
h( x )
H
x
x  x0
 h1 / x1
 h / x
  2 1
 h / x
 m 1
 h1 / xn 
 h2 / xn 





hm / x2  hm / xn 
h1 / x2
h2 / x2
n
m
http://gpspp.sakura.ne.jp
A.4.3.1 最小二乗法
非線形最小二乗法
初期値 x0 の周りで線形化
z  h( x)  ε
 h( x0 )  H ( x  x0 )  ε
最小二乗推定値
ˆx  x0  ( H T WH ) 1 H T W ( z  h( x0 ))
http://gpspp.sakura.ne.jp
A.4.3.2 カルマンフィルタ
逐次推定法
推定値
観測データ
z1 , z 2 ,...
観測更新
xˆ1 , xˆ 2 ,...
時間更新
http://gpspp.sakura.ne.jp
A.4.3.2 カルマンフィルタ
観測更新
K k  Pk () H k ( H k Pk () H k  Rk ) 1
xˆ k ()  xˆ k ()  K k ( z k  h( xˆ k ()))
T
zk
T
xˆ k
Pk ( )  ( I  K k H k ) Pk ()
時間更新
xˆ k 1 ()  xˆ k () 

t k 1
tk
f ( xˆ k (), )d
Pk 1 ()  Φ(t k 1 , t k ) Pk ()Φ(t k 1 , t k )T  Qk
http://gpspp.sakura.ne.jp
A.4.4 相対測位
• 最も一般的な精密測位
• 基本観測量:搬送波位相二重差
• 基準観測点からの相対位置
(基線ベクトル)
• 干渉測位、基線解析
http://gpspp.sakura.ne.jp
A.4.4.1 二重差
LbA
Satellite A
LBa
Satellite B
LBb
LaA
Receiver a
AB
Lab
A
 ( La
Receiver b
B
A
 La )  ( Lb
B
 Lb )
http://gpspp.sakura.ne.jp
A.4.4.1 二重差
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
Lab
AB
AB
AB
AB
AB
  ab  I ab  Tab  nab   ab
AB
AB
AB
  ab  nab   ab (短基
線)
http://gpspp.sakura.ne.jp
A.4.4.2 相対測位の原理
Satellite 3
Satellite 2
Satellite 1
Receiver r
Satellite n
Receiver R
http://gpspp.sakura.ne.jp
A.4.4.2 相対測位の原理
未知パラメータ
12
13
1n T
x  (rr T , N rR
, N rR
,..., N rR
)
x0  (r0T ,0,0,..., 0)T
観測値
zk 
12
13
1n T
( LrR , LrR ,..., LrR )
1
2
  12




/

r



0
r / r0
 rR 
 r
1
3

  13



/

r



r
0
r / r0
rR
hk ( x 0 )  
 Hk  


  
  1 / r   n / r
  1n 
0
r
0
 r
 rR 
0  0 
0   0

   
0 0   

http://gpspp.sakura.ne.jp
A.4.4.2 相対測位の原理
全観測値
z
T
T
T T
 ( z1 , z2 ,..., zm )
h( x0 )  (h1 ( x0 ) , h2 ( x0 ) ,..., hm ( x0 ) )
T
H
T
T
T
 ( H1 , H n ,...,
T T
T T
Hm )
最小二乗推定値
ˆx  x0  ( H T WH ) 1 H T W ( z  h( x0 ))
http://gpspp.sakura.ne.jp
A.4.4.3 整数不定性決定
最小二乗推定値
12 ˆ 13
1n T
ˆ
ˆx  (rˆr T , Nˆ rR
, N rR ,..., N rR )
• 搬送波位相バイアス:実数解
→FLOAT解
• 搬送波位相バイアス:
整数条件を利用した整数化
→FIX解
http://gpspp.sakura.ne.jp
A.4.4.4 相対測位の応用
• RTK(リアルタイムキネマティック)
• 仮想基準点
http://gpspp.sakura.ne.jp
A.4.5 単独測位
軌道/時計
放送暦
~3m/~10 ns
電離層
補正モデル
1~10 m
雑音+マルチパス
コード
0.1~1 m
対流圏
補正モデル
~0.3 m
測位精度
数m~数10m
http://gpspp.sakura.ne.jp
A.4.5 精密単独測位(PPP)
軌道/時計
精密暦
~3cm/~0.1ns
電離層
2周波線形結合
~0.5 cm?
雑音+マルチパス
搬送波位相
0.3~1 cm
対流圏
モデル+推定
~0.5 cm
測位精度
1cm~数cm
http://gpspp.sakura.ne.jp
A.4.5.1 精密単独測位の原理
Satellite 3
Satellite 2
Satellite n
Satellite 1
L2r
L3r
Lnr
L1r
Receiver r
http://gpspp.sakura.ne.jp
A.4.5.1 精密単独測位の原理
未知パラメータ
x
1
2
 (rr , dt, ZW D, N LC r , N LC r ,...,
T
n T
N LC r )
初期推定値・共分散行列
xˆ1 ()  (r0T , dt0 , ZW D0 , N LC 1r 0 , N LC 2r 0 ,..., N LC nr0 )T
P1 ()  diag( x ,  y ,  z ,  dt ,  ZWD ,  N ,...,  N )
2
2
2
2
2
2
2
http://gpspp.sakura.ne.jp
A.4.5.1 精密単独測位の原理
観測更新(1)
1
2
zk  ( LCr , LCr ,...,
n T
LCr )
  1r  c(dtˆ  dT 1 )  Tr1  Nˆ LC 1 
r 

  r2  c(dtˆ  dT 2 )  Tr2  Nˆ LC 2r 
h( xˆ k ())  




  n  c(dtˆ  dT n )  T n  Nˆ n 
r
LC r 
 r
http://gpspp.sakura.ne.jp
A.4.5.1 精密単独測位の原理
観測更新(2)
  1r / rˆr

  r2 / rˆr
Hk  


  n / rˆ
r
 r
c M wet ( Elr1 ) 1 0  0 
c M wet ( Elr2 ) 0 1  0 



   
c M wet ( Elrn ) 0 0  1 
2
2
2
Rk  diag(1 , 2 ,...,  n )
推定値
xˆ k ()  xˆ k ()  K k ( z k  h( xˆ k ()))
http://gpspp.sakura.ne.jp
A.4.5.1 精密単独測位の原理
時間更新
ˆxk 1 ()  (rˆr T , dt0 , ZWˆD, Nˆ LC 1r , Nˆ LC 2r ,..., Nˆ LC nr )T
Pk 1 ()  Φ(t k 1 , t k ) Pk ()Φ(t k 1 , t k )T  QK
Φ(t k 1 , t k )  diag(1,1,1,0,1,1,1,...,1)
Qk  diag(0,0,0,  dt 2 ,  ZWD 2 t ,0,0,..., 0)
http://gpspp.sakura.ne.jp
Up (m)
North (m)
East (m)
A.4.5.1 精密単独測位の例
1
TSKB 2004/10/3 0:00-3:00
0
-1
1
0
-1
1
0
-1
10/3 0:00
10/3 1:00
10/3 2:00
10/3 3:00
http://gpspp.sakura.ne.jp
A.4.5.2 精密単独測位の応用
• キネマティックPPP
• GDGPS(Global Differential GPS)
• 精密時刻同期
• GPS可降水量
http://gpspp.sakura.ne.jp
A.4.6 精密解析の実際
• 解析データ
• 解析データ編集
• 単独測位による概略値計算
• 受信機時計飛びの検出・修正
• サイクルスリップの検出・修正
• パラメータ推定
• アウトライアの検出・除去
http://gpspp.sakura.ne.jp
まとめ
• 精密測位のモデル
幾何学距離、電離層遅延、対流圏遅延、
アンテナ位相中心 ...
• パラメータ推定の手法
• 相対測位
• 精密単独測位(PPP)
• 精密解析の実際
http://gpspp.sakura.ne.jp