ノンパラメトリック分布表現を用いた位置尤度場周辺化による RTK-GNSSの整数アンビギュイティ推定 高野瀬 碧輝，小出 健司，大石 修士，横塚 将志 産業技術総合研究所 (AIST) 

2 国立研究開発法人 産業技術総合研究所 提案手法: サンプリングベースの整数アンビギュイティ推定 RTK-GNSSのFix解推定の新規手法 GPUを利用した広範囲かつ緻密なサンプリングによる整数アンビギュイティ推定手法 キーポイント ✓ サンプリング法を用いた整数アンビギュイティのノンパラメトリック分布推定 ➢ 整数アンビギュイティの不確かさを確率分布によって表現 ➢ 正規分布を仮定しない → マルチパス等の影響を受けても多峰性分布として表現 ✓ GNSS航法演算にGPU実装を導入する新たな試み ➢ GPU並列処理性能を活用した整数アンビギュイティのサンプリング効率の向上 ➢ 提案手法では2700万 (= 3003) サンプルの処理を実行 ( ≒ Fix解) 

3 国立研究開発法人 産業技術総合研究所 提案手法の紹介動画 位置推定結果 サンプリング空間上での尤度分布 Err 整数アンビギュイティの事後確率分布 

4 国立研究開発法人 産業技術総合研究所 提案手法の紹介動画 位置推定結果 サンプリング空間上での尤度分布 Err 整数アンビギュイティの事後確率分布 

5 国立研究開発法人 産業技術総合研究所 提案手法の紹介動画 位置推定結果 サンプリング空間上での尤度分布 Err PAUSE 尤もらしい整数アンビギュイティのみを用いて Fix解を推定する 整数アンビギュイティの事後確率分布 

6 RTK-GNSSの基礎と関連研究 

7 国立研究開発法人 産業技術総合研究所 RTK-GNSSの基礎 RTK-GNSS • 搬送波位相を用いた精密測位の一種 • 波数情報である整数アンビギュイティを求める必要あり • 整数アンビギュイティが求まったら基本的には変動しない →定常的な状態更新が可能 𝜆 𝑁 𝑧 波数 (整数アンビギュイティ) 位相 波長 … 𝑟 𝑥 幾何学距離 状態推定問題 𝑧𝑡 𝑖 = 1 𝜆 𝑟 𝒙𝑡 + 𝑁𝑡 𝑖 状態変数: 𝒙𝑡 ∈ ℝ3 位置 状態変数: 𝑵𝑡 ∈ ℤ𝑘 整数アンビギュイティ 観測変数: 𝒛𝑡 ∈ ℝ𝑘 搬送波位相 𝜆: 波長 𝑟(𝒙𝑡 ): 幾何学距離 整数アンビギュイティと測位ステータス • 整数アンビギュイティが整数解として求まる → Fix解 • 整数アンビギュイティが実数解として求まる → Float解 観測方程式: Fix解を求めるには整数最適化問題を解く必要がある 

8 国立研究開発法人 産業技術総合研究所 関連研究①：LAMBDA法とRatio-test LAMBDA法 (整数最小二乗法) 1. アンビギュイティを実数解として求める (෡ 𝑵𝑡 ∈ ℝ𝑘) 2. 目的関数に従い最適解 ෩ 𝑵𝑡 を探索的に求める ෩ 𝑵𝑡 = arg min 𝑵𝑡 ෡ 𝑵𝑡 − 𝑵𝑡 𝑸𝑁 2 Ratio-test • 推定された整数アンビギュイティの検定 • 閾値以上でFix解として出力 𝑞𝑡 = ෡ 𝑵𝑡 − ෩ 𝑵′𝑡 𝑸𝑁 2 ෡ 𝑵𝑡 − ෩ 𝑵𝑡 𝑸𝑁 2 = 第２最適解の誤差総和 最適解の誤差総和 オープンスカイな環境では高速にFix解探索が可能 アンビギュイティの実数解 ෡ 𝑵𝒕 に依存した目的関数と探索手法 • 実数解 ෡ 𝑵𝑡 周辺に限定した探索 • マルチパスの影響で実数解 ෡ 𝑵𝑡 は精度劣化を引き起こしやすい → 都市部ではFix率が低下する • 局所最適解に陥り，Ratio-testに受かってしまう → ミスFix解の発生 推定した整数アンビギュイティを継続して使用するか，初期化するかは人為的 • サイクルスリップを検出する必要あり 

9 国立研究開発法人 産業技術総合研究所 関連研究②：アンビギュイティフリー手法 AFM/MFAM: (Modified) Ambiguity Function Method • 整数アンビギュイティを解かない手法 • 位置仮説 𝒙𝑡 を用いた残差計算が可能 𝜑 𝑧𝑡 𝑖, 𝒙𝑡 = round 𝑧𝑡 𝑖 − 𝑟(𝒙𝑡 ) 𝜆 − 𝑧𝑡 𝑖 − 𝑟(𝒙𝑡 ) 𝜆 𝜑 𝑧𝑡 𝑖, 𝒙𝑡 : Ambiguity Function Value (AFV) 𝑧𝑡 𝑖:搬送波位相，𝜆:波長，𝑟 𝒙𝑡 :幾何学距離 MU-PF: Multiple Update Particle Filter • AFMの多峰性分布を多層リサンプリングで単峰性分布近似 少数のサンプリングでもAFMの多峰性分布から最適解を抽出可能 マルチパスの影響で尤度を一意に決定できない影響あり → 多重仮説性× サンプリングする位置範囲が限定的 or サンプリング密度が疎 • 真の整数アンビギュイティ周辺にパーティクルが収束しない可能性あり 𝑝 𝒛𝑡 |𝑵𝒕 , 𝒙𝑡 → 𝑝 𝒛𝑡 |𝒙𝑡 サイクルスリップ問題に頑強 正確な分布表現には膨大なサンプル数が必要 波長由来 (L1:0.19m) 

10 提案手法 

11 国立研究開発法人 産業技術総合研究所 提案手法の概要 GPUを利用した広範囲かつ緻密なサンプリングによる整数アンビギュイティ推定手法 アプローチ ✓ サンプリング法を用いたノンパラメトリック分布推定 ✓ ヒストグラムフィルタを用いた確率的な状態更新 課題: 整数アンビギュイティ空間上でのサンプリングは難しい • 真値を含む分布の定義域が不明 • 衛星間の相互関係のモデル化が困難 → サンプリング効率が向上できない • サンプル間の組合せ爆発 → 衛星数k個で，1k, 10k, 20k, …個の組合せパターン 提案手法: 位置尤度場周辺化 • AFMを応用して位置空間上で整数アンビギュイティの確率分布をサンプリングする 膨大なサンプル数が必要 → 3軸で，1003, 2003, 3003, …個 GPUによって計算処理速度を向上させる 𝑝 𝑵𝑡 |𝒛𝑡 ∝ 𝑝 𝒛𝑡 |𝑵𝑡 𝑝 𝑵𝑡 事後分布 尤度分布 事前分布 状態変数: 𝑵𝑡 ∈ ℤ𝑘 整数アンビギュイティ 観測変数: 𝒛𝑡 ∈ ℝ𝑘 搬送波位相 

12 国立研究開発法人 産業技術総合研究所 位置尤度場周辺化による整数アンビギュイティの尤度分布推定 位置尤度場周辺化とは？ 𝑝 𝒛𝑡 |𝑵𝑡 = න 𝑝 𝒛𝑡 |𝑵𝑡 , 𝒙𝑡 𝑝 𝒙𝑡 𝑑𝒙𝑡 尤度分布 AFMによる𝑵𝒕 の変数消去 𝑝 𝒛𝑡 |𝑵𝑡 = න 𝑝 𝒛𝑡 |𝒙𝑡 𝑝 𝒙𝑡 𝑑𝒙𝑡 尤度分布 𝒙𝑡 と 𝒛𝑡 だけで 𝑝 𝒛𝑡 |𝑵𝑡 が計算可能になる！ 状態変数: 𝒙𝑡 ∈ ℝ3 位置 状態変数: 𝑵𝑡 ∈ ℤ𝑘 整数アンビギュイティ 観測変数: 𝒛𝑡 ∈ ℝ𝑘 搬送波位相 提案手法: 位置尤度場周辺化による整数アンビギュイティの尤度分布推定 位置空間の 化 A Mによる の尤度 位置尤度場の 位置尤度場の 周辺化 

13 国立研究開発法人 産業技術総合研究所 位置空間の 化 A Mによる の尤度 位置尤度場の 位置尤度場の 周辺化 任意の位置空間の 化 𝑥 𝑦 𝑧 • 各サンプルは一様分布 • 任意の位置空間を格子状に離散化する 

14 国立研究開発法人 産業技術総合研究所 の尤度 𝑥 𝑦 𝑧 • 各サンプルの尤度分布を計算 𝑝 𝑧𝑡 𝑖|𝒙𝑡 = 1 2𝜋𝜎𝑧 2 exp − 1 2 𝜑 𝑧𝑡 𝑖, 𝒙𝑡 2 𝜎𝑧 2 位置空間の 化 A Mによる の尤度 位置尤度場の 位置尤度場の 周辺化 A Mを利用した尤度関数 𝜑 𝑧𝑡 𝑖, 𝒙𝑡 : Ambiguity Function Value (AFV) 𝑧𝑡 𝑖:搬送波位相，𝒙𝑡 :位置，𝜎𝑧 2: ノイズパラメータ 

15 国立研究開発法人 産業技術総合研究所 の尤度 𝑥 𝑦 𝑧 位置空間の 化 A Mによる の尤度 位置尤度場の 位置尤度場の 周辺化 • 各サンプルの尤度分布を計算 

16 国立研究開発法人 産業技術総合研究所 位置尤度場の • 全衛星の尤度分布の総積を計算する 位置空間の 化 A Mによる の尤度 位置尤度場の 位置尤度場の 周辺化 𝑝 𝒛𝑡 |𝒙𝑡 = ෑ 𝑖 𝑝 𝑧𝑡 𝑖|𝒙𝑡 →位置尤度場 𝑝 𝑧𝑡 𝑖|𝒙𝑡 = 1 2𝜋𝜎𝑧 2 exp − 1 2 𝜑 𝑧𝑡 𝑖, 𝒙𝑡 2 𝜎𝑧 2 

17 国立研究開発法人 産業技術総合研究所 位置尤度場の周辺化 位置空間の 化 A Mによる の尤度 位置尤度場の 位置尤度場の 周辺化 • 位置尤度場を周辺化 𝜑 𝑧𝑡 𝑖, 𝒙𝑡 = round 𝑧𝑡 𝑖 − 𝑟(𝒙𝑡 ) 𝜆 − 𝑧𝑡 𝑖 − 𝑟(𝒙𝑡 ) 𝜆 各位置サンプルに対する 𝑁𝑡 に相当 整数アンビギュイティ𝑁𝑡 との対応付けは？ 尤度を𝑵𝒕 のビンの重みとして加算 

18 国立研究開発法人 産業技術総合研究所 位置尤度場の周辺化 位置空間の 化 A Mによる の尤度 位置尤度場の 位置尤度場の 周辺化 • 位置尤度場を周辺化 𝜑 𝑧𝑡 𝑖, 𝒙𝑡 = round 𝑧𝑡 𝑖 − 𝑟(𝒙𝑡 ) 𝜆 − 𝑧𝑡 𝑖 − 𝑟(𝒙𝑡 ) 𝜆 各位置サンプルに対する 𝑁𝑡 に相当 整数アンビギュイティ𝑁𝑡 との対応付けは？ 尤度を𝑵𝒕 のビンの重みとして加算 

19 国立研究開発法人 産業技術総合研究所 位置尤度場の周辺化 位置空間の 化 A Mによる の尤度 位置尤度場の 位置尤度場の 周辺化 • 位置尤度場を周辺化 整数アンビギュイティ 確立密度 

20 国立研究開発法人 産業技術総合研究所 位置尤度場の周辺化 位置空間の 化 A Mによる の尤度 位置尤度場の 位置尤度場の 周辺化 • 位置尤度場を周辺化 整数アンビギュイティ 確立密度 整数アンビギュイティ 確立密度 整数アンビギュイティ 確立密度 Satellite 3 Satellite 2 Satellite 1 

21 国立研究開発法人 産業技術総合研究所 ヒストグラムフィルタを用いた状態更新 𝑝 𝑵𝑡 |𝒛𝑡 ∝ 𝑝 𝒛𝑡 |𝑵𝑡 𝑝 𝑵𝑡 事後分布 尤度分布 事前分布 整数アンビギュイティの事後分布推定 = න 𝑝 𝒛𝑡 |𝒙𝑡 𝑝 𝒙𝑡 𝑑𝒙𝑡 න 𝑝 𝑵𝑡 |𝑵𝑡−1 𝑝 𝑵𝑡−1 |𝒛𝑡−1 𝑑𝑵𝑡−1 尤度分布 事前分布 (定常状態モデル) 事前分布と尤度分布との関連性 ✓ 形状が一致 → 事後確率が収束する ✓ 形状が不一致 ＝ 尤度分布がマルチパス/サイクルスリップの影響で変化する → 不確実性が保持される ＊事前分布の初期値は一様分布 × ∝ × ∝ 事前分布 尤度分布 事後分布 意図的にノイズ付与 

22 国立研究開発法人 産業技術総合研究所 整数アンビギュイティの代表値抽出と位置推定 代表値抽出: 最大事後確率による抽出 位置推定: ガウスニュートン法 ✓ 本研究は事後確率の尤もらしさを閾値によって判定する ✓ 今後はノンパラメトリック検定によるに移行予定 (e.g.コルモゴロフ–スミルノフ検定) ✓ 簡単な外れ値対策のため、Huberカーネル 𝜌 によるロバスト最適化を適用する ෩ 𝑵𝑡 = arg max 𝑵𝑡 𝑝 𝑵𝑡 |𝒛𝑡 ෥ 𝒙𝑡 = arg min 𝒙𝑡 ෍ 𝑖=1 𝑘 𝜌 𝑧𝑡 𝑖 − 1 𝜆 𝑟 𝑥𝑡 + ෩ 𝑁𝑡 𝑖 2 

23 実験 静止体環境/移動体環境 

24 国立研究開発法人 産業技術総合研究所 提案手法のパラメータ (整数アンビギュイティの事後分布推定) 分割数 分解能 Res. ノイズパラメータ Noise 位置空間のサイズとサンプル数が決まる e.g. 分解能 0.02m, 分割数 100 → 各軸 2m, サンプル数 1003 e.g. 分解能 0.10m, 分割数 300 → 各軸 30m, サンプル数 3003 位置空間の各軸 • 分解能 Res. • 分割数 ノイズパラメータ Noise 𝜎𝑧 2 • AFMによる尤度関数の重みパラメータ 𝑝 𝑧𝑡 𝑖|𝒙𝑡 = 1 2𝜋𝜎𝑧 2 exp − 1 2 𝜑 𝑧𝑡 𝑖, 𝒙𝑡 2 𝜎𝑧 2 

25 国立研究開発法人 産業技術総合研究所 静止体環境の実験設定 実験データ • 基準局データ(つくば1, つくば3) 実験1: パラメータ毎の位置推定性能評価 • 位置尤度場の中心は真値 • 各軸 100分割 = 1003サンプル • 分解能とノイズ の設定は以下の表通り 実験2: CPU実装 vs GPU実装 の処理時間評価 • 位置尤度場の中心は真値 • 各軸100分割，200分割，300分割 • 分解能: 0.02m，ノイズ: 0.001m2 AFMによる尤度関数 分割数 分解能 Res. 提案手法のパラメータ概要 ノイズ Noise 使用計算機 • CPU: Intel Core i9-14900KF • GPU: NVIDIA GeForce RTX 4080 

26 国立研究開発法人 産業技術総合研究所 静止体環境の実験結果 分解能(高) ✓ 分解能が低いと位置推定精度が劣化する ➢ 尤度分布の正確性が低下する ➢ 事後確率が一意に収束しなくなる ✓ GPUによる処理時間の高速化 分解能(低) 実験2: 処理時間結果 実験1: 位置推定結果 整数アンビギュイティ 確立密度 

27 国立研究開発法人 産業技術総合研究所 移動体環境の実験設定 実験データ • OpenData (Suburban, Urban 1) • 高精度測位チャレンジ2023 (Urban2, Urban3) 使用計算機 • CPU: Intel Core i9-14900KF • GPU: NVIDIA GeForce RTX 4080 実験1: LAMBDA法との位置推定性能比較評価 • ベースライン: RTKLIB, RTKLIB demo5 • 位置尤度場の中心はFloat解(位置) • 分解能: 0.02m，ノイズ: 0.001m2 • 各軸 300分割 = 3003サンプル 実験2: 整数アンビギュイティ事後分布の定性評価 • 環境に応じた事後分布の調査 • 実験1 Urban 3 の結果を使用 Float解を中心とした 各軸6mの位置空間を密にサンプリング AFMによる尤度関数 分解能 Res. 提案手法のパラメータ概要 ノイズ Noise 

28 国立研究開発法人 産業技術総合研究所 移動体環境の実験結果 Float解の誤差平均 誤差0.5m以内のCDF ✓ Float解の誤差が小さい → LAMBDA法ベースの手法と同程度な推定結果 ✓ Float解の誤差が大きい → 提案手法が優位 サンプリング空間の広域化によって真値を含む確率を向上させた 実験1: 位置推定結果 

29 国立研究開発法人 産業技術総合研究所 移動体環境の実験結果 Float解の誤差平均 誤差0.5m以内のCDF ✓ Float解の位置精度が悪く、サンプル空間外に真値を多く含む場合 →提案手法ではFix解を推定することは困難 →サンプル空間の拡大化などの対策が必要 (今後の課題) 実験1: 位置推定結果 

30 国立研究開発法人 産業技術総合研究所 移動体環境の実験結果 Point A 誤差: 0.19m Point B 誤差: 0.25m Point C 誤差: 1.68m Point D 誤差: 1.96m 誤差 [m] 0.0 3.0 G10 Point A 誤差: 0.19m G25 Point A 誤差: 0.19m G29 Point A 誤差: 0.19m Point B 誤差: 0.25m G25 Point B 誤差: 0.25m G29 Point B 誤差: 0.25m G10 Point C 誤差: 1.68m G25 Point C 誤差: 1.68m G29 Point C 誤差: 1.68m G10 Point D 誤差: 1.96m G25 Point D 誤差: 1.96m G29 Point D 誤差: 1.96m 整数アンビギュイティ 確率密度 ✓ 分布が収束する → 位置誤差が小さい ✓ 分布が多峰性である → 整数アンビギュイティが一意に決まらず，位置誤差が大きい Fix解を推定しにくい環境であることの裏付け 実験2: 整数アンビギュイティ事後分布の定性評価 

31 今後の課題 

32 国立研究開発法人 産業技術総合研究所 ix解のオフライン解析 提案手法の整数アンビギュイティ推定は単純なベイズフィルタリング ✓ 容易にバックワードスムージングの定式化が可能 ✓ より尤もらしいFix解推定への貢献に期待 Fix解が取得できなかった環境での定性評価の可視化 ✓ 「なぜFixしないのか」の説明を可能に！ →Float解周辺のサンプリングで 整数アンビギュイティの尤度分布が分かる ✓ Fix解の信頼性の評価に活用？ ✓ 分布を活用したマルチパス除去？ 

33 国立研究開発法人 産業技術総合研究所 サンプリング性能について サンプリングする位置空間の大きさ/密度と位置精度の関係性 位置空間を広範囲かつ緻密に離散化すると？ 位置空間の中心位置精度にゆとりができる ゆっくり移動体では高頻度に位置空間の移動更新が不要 計算コストが指数関数的に増加する 局所最適解に陥る可能性が高くなる 多様なセンサとの複合が可能になる ✓ 提案手法は位置空間上のサンプリングでFix解を推定できる！ ✓ 位置空間は任意なので独立していてもよい！ ✓ 各サンプリング空間を縮小させサンプリング性能の向上 (スマホ等の運動予測が困難なアプリに応用可) x y x y Grand truth GNSS Solution Sampling area 

34 まとめ 

35 国立研究開発法人 産業技術総合研究所 まとめ 整数アンビギュイティを推定する恩恵 Fix解が得られる 整数アンビギュイティの定常特性を利用したい →正しい推定/状態更新で、高精度位置を維持したスムージング 整数アンビギュイティのノンパラメトリック分布表現 ＆ ベイズ更新の確立 整数アンビギュイティの確率的な状態更新 ノンパラメトリック分布表現によってマルチパスやサイクルスリップへの頑強化 位置空間上でのサンプリングで整数アンビギュイティの尤度分布推定 整数アンビギュイティの直接的なサンプリングが困難な課題の解決 GPU並列処理を活用してサンプリング処理能力の強化