[Peer Reviewed]
Proceedings of the 4th Conference of the Society for the Advancement of Positioning Technologies (SAPT), Osaka, Japan, Sep. 9, 2022

サイクルスリップの検出と修復による PPP の性能向上
Improvement of PPP performance by detecting and repairing cycle slips
山内義洋，辻井利昭

大阪公立大学

Yoshihiro Yamauchi，Tsujii Toshiaki
Osaka Metropolitan University
Abstract
A major issue with PPP is convergence time. When the GNSS signal is interrupted or attenuated,a cycle slip might occur, which results in the initialization of the Ambiguity estimates. As a result,additional time to re-converge the estimates is necessary each time a cycle slip occurs. The purposeof this research is to prevent initialization of the estimates by detecting the cycle slip and removing
it from the observed value. A simulation study using adopted algorithm has shown that the cycleslip is detected and repaired correctly.
１．研究背景・目的
近年，自動運転が開発，導入されつつある中で，その完全な実現には数𝑐𝑚～数 10𝑐𝑚といった高精度の測位が要求される．そこで精密単独測位方式 (図 1 PPP：Precise Point Positioning) を利用する．相対測位では近傍基準局との観測値の差を取ることで高精度な測位を可能としている．一方で PPP は搬送波位相を主に利用し精密歴は外部データから取り込むことにより，近傍の基準局を必要とせずにあらゆる場所において誤差数 cm～数 10cm の測位が可能である．

図1

PPP と相対測位の概念図

しかしながら衛星信号を遮断・反射する都市部や山中では頻繁にサイクルスリップが生じやすい．サイクルスリップとは信号の攪乱により，通常一定のアンビギュイティ𝑁が不連続に変化してしまう事である．発生要因としては，遮蔽環境による電力損失及び低い信号対雑音比や衛星発振器の誤作動，電離層の擾乱が挙げられる．図 2 はアンビギュイティが不連続に変化する様子を示している．PPP では測位
収束に 10～20 分要するが，サイクルスリップが生じると推定値をリセットするため，再び収束時間が必要となるという課題がある．そこで不連続な変化量∆𝑁を検出し搬送波観測量から除去する事により，推定値の初期化が引き起こす測位精度の劣化を防ぐ事を目的とした．

図2

第4回測位技術振興会研究発表講演会
講演論文集（2022年9月9日，大阪）

サイクルスリップ発生図-40-

２．精密単独測位(PPP)の概要
まずは観測量によって得られる方程式及びそれらの時間差分方程式を式(1)~ (4)に示す．
𝑘𝛷𝑖𝑘 = 𝜌𝑖𝑘 + 𝑐𝑑𝑡𝑟𝑘 − 𝑐𝑑𝑡𝑠𝑘 − 𝐼𝑖𝑘 + 𝑇 𝑘 + 𝜆𝑖 (𝑁𝑖𝑘 + 𝑙𝑖𝑘 ) + 𝜀∅𝑖(1)𝑘
𝑃𝑖𝑘 = 𝜌𝑖𝑘 + 𝑐𝑑𝑡𝑟𝑘 − 𝑐𝑑𝑡𝑠𝑘 + 𝐼𝑖𝑘 + 𝑇 𝑘 + 𝜀𝑃𝑖(2)𝑘∆ 𝛷𝑖𝑘 = ∆𝜌́ 𝑖𝑘 − ∆𝐼𝑖𝑘 + 𝜆𝑖 ∆𝑁𝑖𝑘 + ∆𝜀∅𝑖

(3)𝑘∆𝑃𝑖𝑘 = ∆𝜌́ 𝑖𝑘 + ∆𝐼𝑖𝑘 + ∆𝜀𝑃𝑖

(4)

𝑘：衛星番号
𝑖：L1,L2,L5 信号(𝑖 = 1,2,3)∆：時間差分𝛷：搬送波位相[𝑚]𝑃：擬似距離 [𝑚]𝜌：衛星と受信機間の幾何学的距離[𝑚]𝜌́ ：衛星と受信機間の幾何学的距離+時計誤差+対流圏遅延を含む換算量[𝑚]𝑐𝑑𝑡𝑟 ：受信機時計誤差[𝑚]𝑐𝑑𝑡𝑠 ：衛星時計誤差[𝑚]𝑇：対流圏遅延量[𝑚]𝐼：電離層遅延量[𝑚]𝑁：アンビギュイティ[𝑐𝑦𝑐𝑙𝑒]
𝑙：衛星初期位相バイアス[𝑐𝑦𝑐𝑙𝑒]
𝜆：波長[𝑚]𝜀：マルチパス含む観測ノイズ[𝑚]PPP では式(1),(2)において，精密歴を利用し𝑐𝑑𝑡𝑟 , 𝑐𝑑𝑡𝑠 及び 𝑙を補正し，𝐼は電離層フリー線形結合により排除し，𝑇はモデルにより推定を行う．𝜌については単独測位解を初期値とし，𝑁は電離層フリー線形結合により得られる値を初期値として，カルマンフィルタによる観測更新・時間更新により推定し，測位を求めている．本研究では 4 章で述べるサイクルスリップの処理を行った後に PPP を利用している．
３．サイクルスリップによる測位劣化
サイクルスリップが生じた際に測位精度に及ぼす影響について検証を行った．今回の検証ではサイクルスリップを人為的に発生させたデータを利用する．そのためサイクルスリップが生じにくいオープン
スカイ環境で静止データを取得し，特定の時刻の搬送波観測値を整数サイクル増減させることでサイクルスリップを発生させた．実験は 2021 年 6 月 28 日に大阪公立大学中百舌鳥キャンパス B6 棟屋上で行い，図 3 及び図 4 ではその実験場所とその時の衛星配置を示しており，赤色の衛星軌道は L1 帯と L2 帯信号を発信する衛星を，青色の衛星軌道はそれらに加えて L5 帯信号を発信する衛星を表している．

-41-

図 3 実験場所(大阪公立大学 B6 棟屋上)

図4

衛星配置

なお、解析ソフトとして、GNSS 解析ソフトウェア RTKLIB [2]を使用した．PPP 及び相対測位に関して，RTKLIB の主要な設定はそれぞれ表 4,5 の通りである．この設定はこの後のすべてのデータ解析で同一である．また PPP で必要な高精度クロック・軌道及び衛星位相バイアスの補正値は CNES が提供するリアルタイム・プロダクトを使用した．[3]
相対測位は収束がほぼ瞬時であり常時 cm 級の精度を保つことから真値として定め，基準局は国土地理院の堺電子基準点（大阪府）とし，そのオープンデータ[4]を用いた．
表 1 PPP における RTKLIB の設定項目
表 2 相対測位における RTKLIB の設定項目
設定項目

設定値

設定項目

設定値

解析ソフトウェア

RTKLIBver2.4.2master

解析ソフトウェア

RTKLIBver2.4.2master

測位方式

後処理 PPP(kinematic)

測位方式

後処理 static

利用衛星

GPS

利用衛星

GPS,QZSS

利用周波数帯

L1,L2,L5

利用周波数帯

L1,L2,L5

衛星精密歴ファイル

cnt.sp3,cnt.clk(CNES)

衛星仰角マスク15°

衛星位相バイアスファイル

cnt.bia(CNES)

基準点観測データ

堺(電子基準点)

衛星仰角マスク

15°

電離層補正

Broadcast

電離層補正

Iono-free LC

対流圏補正

Saastamoinen

対流圏補正

Estimate ZTD

実験は GPST05:15~07:15 の約 2 時間行い，G09 衛星に対し GPST06:30 にサイクルスリップを発生させた．実験開始後約１時間に設定したのは PPP の収束時間を十分に設けるためであり，サイクルスリップ量は(∆𝑁1 , ∆𝑁2 , ∆𝑁3 ) = (1,1, 1)とした．この表記は L1,L2,L5 信号の搬送波観測量を全て 1cycle 増加させたことを表している．またマルチパス環境を想定し，観測衛星数による以下 2 ケース(i),(ii)の比較を行った． それぞれの衛星配置は図 5,7 に示し，真値に対する PPP の測位精度は図 6,8 に示している．

-42-

(i)全観測衛星(約 10 機)を利用した場合

図5衛星配置

図6

PPP の測位精度

図8

PPP の測位精度(ii)高仰角衛星を 5 機利用した場合

図7衛星配置

利用衛星数が約 10 機の場合に対し 5 機と少ない場合，サイクルスリップが生じた時刻で測位結果が不連続に跳んでしまうことが分かる．この理由はサイクルスリップを検知した際に RTKLIB ではアンビギュイティの推定値の初期化を行うからであり，利用衛星数が少ないほど測位に与える影響が相対的に大きくなってしまうためである．したがって可視衛星が限定されるマルチパス環境において，サイクルスリップが生じやすいだけでなく，推定値の初期化による測位精度の大幅な劣化を引き起こす可能性がある．
そのため PPP では測位の再収束と精度劣化を防ぐため，サイクルスリップを検出するだけでなく除去を行う必要があると考えられる．次章ではそのサイクルスリップの処理について述べる．
４．サイクルスリップ処理アルゴリズム
サイクルスリップの処理を①検出②推定③修復の三段階で行い，搬送波位相観測量を修正した後に 2章で述べた精密単独測位を行うこととする．通常，𝑁は一定値であり∆𝑁 = 0であるが，サイクルスリップが生じた場合，整数値を示す．
① サイクルスリップの検出
文献[1]に基づき，式(5)のサイクルスリップ検出式を導入した．

-43-

|𝑤1 ∆∅1 + 𝑤2 ∆∅2 + 𝑤3 ∆∅3 |222√2√𝑤1 2 ∙ 𝜆1 + 𝑤2 2 ∙ 𝜆2 + 𝑤3 2 ∙ 𝜆3

≥ 𝑓 ∙ 𝜎𝐿1𝑐𝑦𝑐𝑙𝑒(5)

𝜎𝐿1𝑐𝑦𝑐𝑙𝑒 ：L1 観測ノイズ(0.05𝑐𝑦𝑐𝑙𝑒)
𝑓：観測周波数(= 1,2,3)
また𝑤1 + 𝑤2 + 𝑤3 = 0を満たし，𝑤1 , 𝑤2 , 𝑤3 は整数である．
左辺分子の式は 3 周波信号を利用したジオメトリフリー線形結合を示している．∆𝜌の係数を 0 にでき，周波数に依存しない項は排除できるため，電離層遅延量の時間差分∆𝐼と∆𝑁項のみで構成される．そのため∆𝑁を的確に検知するためには∆𝐼を小さくする必要がある．表 3 は文献[1]を参照した(𝑤1 , 𝑤2 , 𝑤3 )の組に対する∆𝐼の大きさについて示している．∆𝐼𝑖 の影響を小さくする組(𝑤1 , 𝑤2 , 𝑤3 ) = (−1, −1, 2), (−1, 4, −3)と定めることで, ∆𝐼を最小にでき∆𝑁項が主に左辺を支配するパラメータとみなすことができる．これらの組み合わせのどちらか一方が不等式(5)を満たせば検出とみなした．(図 9) 通常𝜎𝐿1𝑐𝑦𝑐𝑙𝑒 は0.01𝑐𝑦𝑐𝑙𝑒程だが，マルチパスノイズを含めるため0.05𝑐𝑦𝑐𝑙𝑒と設定している．
表 3 サイクルスリップ検出の為の 2 段階ジオメトリフリーの組み合わせ
第 1 段階の組(𝑤1 , 𝑤2 , 𝑤3 )第 2 段階の組(𝑤1 , 𝑤2 , 𝑤3 )(−1, −1, 2)(−1, 4, −3)∆𝐼の合計量1.7414∆𝐼(−1,2, −1)(−1, −3, 4)

1.9639∆𝐼(−1, −1, 2)(−1, 3, −2)

1.8295∆𝐼(−1,2, −1)(−1, −2, −3)

2.0305∆𝐼

図9

検出式の簡略図② サイクルスリップの推定式(3)を変形し，計 2 回時間差分を取った以下の式(6)~(9)を導入し，∆𝑁の推定を行った．高レートのデータを活用する事で観測値及び誤差推定値の時間変化量が一定とみなすことができ，差分により排除できることを利用した方法で，また本稿ではこの推定方法を 2 回差分方式と呼んでいる．
∆∆𝛷 𝑖 = ∆∆𝜌́ − ∆∆𝐼𝑖 + 𝜆𝑖 ∆∆𝑁𝑖 + ∆∆𝜀 𝛷𝑖(6)3 エポックを利用し式を展開すると，式(7)が得られる．また電離層遅延量と観測ノイズの 2 回差分項は 0と近似した．
𝑡𝑡𝑡𝑡(∅𝑖 0 − ∅𝑖 −1 ) − (∅𝑖 −1 − ∅𝑖 −2 )𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡= (𝜌́ 𝑖 0 − 𝜌́ 𝑖 −1 ) − (𝜌́ 𝑖 −1 − 𝜌́ 𝑖 −2 ) + 𝜆𝑖 (𝑁𝑖 0 − 𝑁𝑖 −1 ) − 𝜆𝑖 (𝑁𝑖 −1 − 𝑁𝑖 −2 )𝑡0 , 𝑡−1 , 𝑡−2 ：利用エポック-44-(7)ここでサイクルスリップが起きていない場合では，通常𝑁は一定であるため，𝑁項の 2 回差分は 0 となる．また図 10 のように搬送波観測量の時間変化は高レートによる直線近似により一定とみなせるため，左辺は 0 と近似でき，残された右辺の幾何学的距離項も 0 となる．

図 10

搬送波位相観測量の時間変化
𝑡𝑡𝑡𝑡

次にサイクルスリップが時刻𝑡0 で発生したと仮定すると，∆𝑁𝑖 = 𝑁𝑖 0 − 𝑁𝑖 −1 , 𝑁𝑖 −1 − 𝑁𝑖 −2 = 0と表せ，式(8)が得られる．
𝑡𝑡𝑡𝑡(∅𝑖 0 − ∅𝑖 −1 ) − (∅𝑖 −1 − ∅𝑖 −2 ) = 𝜆𝑖 ∆𝑁𝑖(8)したがって，高レートのデータを利用する事で∆𝑁𝑖 を以下の式で導出できる．
𝑡𝑡𝑡(9)∆𝑁𝑖 = (∅𝑖 0 − 2∅𝑖 −1 − ∅𝑖 −2 )/𝜆𝑖

今回の解析では 2 回差分の観測値及び誤差推定値を 0 と近似できる 20Hz のデータを利用している．
③ 搬送波位相観測値の修復求めた∆𝑁1 , ∆𝑁2 , ∆𝑁3 を用いて，式(10)~(13)により各搬送波位相観測値の修復を行う．そして得られた観測値がサイクルスリップを排除できているか不等式(13)により再評価する．また不等式(13)は不等式(5)に対し不等号の向きが反対である．
𝑟𝑒𝑝𝑎𝑖𝑟= ∆∅1 − 𝜆1 𝛥𝑁1(10)𝑟𝑒𝑝𝑎𝑖𝑟= ∆∅2 − 𝜆2 𝛥𝑁2(11)𝑟𝑒𝑝𝑎𝑖𝑟= ∆∅3 − 𝜆3 𝛥𝑁3(12)∆∅1∆∅2∆∅3𝑟𝑒𝑝𝑎𝑖𝑟|𝑤1 ∆∅1𝑟𝑒𝑝𝑎𝑖𝑟+ 𝑤2 ∆∅2𝑟𝑒𝑝𝑎𝑖𝑟+ 𝑤3 ∆∅3222√2√𝑤1 2 ∙ 𝜆1 + 𝑤2 2 ∙ 𝜆2 + 𝑤3 2 ∙ 𝜆3|≤ 𝑓 ∙ 𝜎𝐿1𝑐𝑦𝑐𝑙𝑒(13)したがって不等式(13)を満たせば，図 11 のように閾値を下回り，∆𝑁を適切に除外できたとみなせ，アンビギュイティを一定に保つ事ができる．

図 11

修復式の簡略図-45-

５．アルゴリズム性能評価
4 章の考察より，サイクルスリップの影響を大きく受ける可視衛星数が少ないマルチパス環境のデータを解析に用いた．実験場所は大阪公立大学中百舌鳥キャンパス B2 棟 B3 棟間(図 12)で，2022 年 9 月 6 日に衛星信号を約 10 分間 20Hz で取得した．
５.１ シミュレーション
今回は人為的に発生させたサイクルスリップの推定ができるのかを実証するため，衛星配置から建物に遮蔽されることなく，受信し続けることができた衛星のみを利用する．その時のスカイプロットを図 13に示す．このデータを用いてサイクルスリップ処理アルゴリズムの検証を行った．そしてサイクルスリッ
プ処理後は同一の設定で PPP を算出し，処理による測位精度を比較する事で性能の評価を行った．

図 12 実験場所(大阪公立大学 B6 棟屋上)

図 13

衛星配置

シミュレーションとして，取得した搬送波位相観測値に人為的な∆𝑁を以下 3 パターンに区別し加えた．
(a)微小な∆𝑁 (∆𝑁 ≤ 1)(b)3 周波のそれぞれが等しい∆𝑁(c)比較的大きな∆𝑁 (∆𝑁 > 1)
そこで先述したアルゴリズムにより，∆𝑁を検出・推定し，その性能を確かめた．
５.２ ∆𝑵推定結果
推定値∆𝑁と人為的に発生させた∆𝑁を表 4 で比較した．
例)1 行目：15𝑒𝑝𝑜𝑐ℎ目の G06 衛星の L1,L2 搬送波位相観測値を 1𝑐𝑦𝑐𝑙𝑒ずつ増加
表 4 各パターンにおける推定値∆𝑁-46-今回の実証実験では意図的に発生させた∆𝑁を検出・推定できていたため，∆𝑁の処理は適切に行えていたと言える．
５.３ 静止体測位
次にサイクルスリップ検出・推定・修復による収束改善について考察する．可視衛星数 4 機の中から最も高仰角である G24 に着目し，GPST07:30 にサイクルスリップを発生させ，サイクルスリップの処理を行った後に測位結果を求めた． 図 14 はサイクルスリップの処理を行わない場合を示し，図 15 は本アルゴリズムを導入した結果を表している．

図 14

PPP の測位結果(∆𝑁未処理)

図 15

PPP の測位結果(∆𝑁処理)

５.４ 結果・考察
アルゴリズムの性能評価の為，様々なパターンのサイクルスリップを発生させ検証を行った結果，正確な値を推定する事ができていた．またサイクルスリップによる測位劣化については図 14,15 より，サイクルスリップを除去したことで推定値の初期化及び測位の再収束を防ぐことができたと言える．本来高精度である PPP の測位誤差が数 m に及んでいる原因としては，観測衛星数が少ないことにより測位収束が完全に終えていないからだと考えられる．
６．結論
3 章ではサイクルスリップがどのような状況で測位に甚大な影響を及ぼすのかを検討し，観測衛星数が極端に減るマルチパス環境下であると考えられた．その理由としてはアンビギュイティ推定値の初期化による影響が相対的に大きくなってしまう事がと考えられる．4 章で述べたサイクルスリップ検出・推定・-47-

除去アルゴリズムの評価を 5 章で行った．その結果，静止点観測データに人為的な∆𝑁を発生させたタイミングで∆𝑁を検知し，正しく推定し，搬送波位相から取り除くことができていると確認された．
７．今後の課題
今回の研究内容は実データに人為的なサイクルスリップを加えた状況におけるサイクルスリップアルゴリズムの性能評価であったため，次の段階としてサイクルスリップが自然発生した実データでサイクルスリップを検知・推定・除去できるか検討する．また現段階では GPS の 3 周波信号を観測できた衛星しか利用していないため，今後は 2 周波信号や他の GNSS に対応させていく．そして最終的な目標としては移動体データに適用し，特にマルチパス環境における測位収束の改善，再初期化の防止を取り組む予定である．
参考文献[1] Zhen Dai,Stefan Knedlik and Otmar Loffeld “Instantaneous Triple-Frequency GPS Cycle-Slip Detection andRepair“ ,International Journal of Navigation and Observation ,Vol 2009 ,Article ID 407231[2] 高須知二『RTKLIB: An Open Source Program Package for GNSS Positioning』
http://www.rtklib.com/(2023 年 1 月 30 日閲覧)
[3] CNES『The PPP-WIZARD project Real-time products』
http://www.ppp-wizard.net/products/REAL_TIME/(2023 年 1 月 30 日閲覧)
[4] 国土地理院 『電子基準点データ提供サービス』 https://terras.gsi.go.jp/(2023 年 1 月 30 日閲覧)

-48-