生のGPSデータを用いた位置計算 (2017)
(telesens.co)- GPSの位置計算は、衛星までの 疑似距離(pseudorange) と軌道暦、受信機時計誤差を同時に解く問題であり、生データをMatlabで処理して実際の位置推定まで実装している
- 計算には緯度・経度よりも WGS 84 ECEF 座標系が適しており、ユーザー基準の方位角・仰角や誤差解析には ENU 局所座標系も併用される
- 衛星位置はGPS Interface Specificationの手順と軌道パラメータから求め、信号伝搬中に地球が回転するため、送信時点の座標を受信時点のECEFフレームに合わせて補正する必要がある
- ユーザー位置と受信機時計バイアスは、少なくとも4個の衛星の補正済み疑似距離に対して 反復最小二乗法 を行って推定し、例の解析では電離層・対流圏遅延を除外している
- u-blox NEO-M8T/6T、RTKLib STRSVR、RTCM 1002/1019、goGPSを用いた静止受信機実験では、位置の標準偏差は East 14.00m、North 39.88m、Up 47.35mで、時計バイアスは 4.27e-7sec/sec で線形にドリフトした
GPS位置計算の基本問題
- GPSの中核的な作業はユーザーの 位置 を計算すること
- 緯度・経度・高度は地表位置の表現にはなじみ深いが、緯度や経度1度の差が常に同じ物理距離を意味するわけではないため、数学的計算には不便
- 経度1度差の距離は赤道で最大となり、極地では0に近づく
- 計算には、単位座標差が一定の物理距離を表す直交座標系が必要
- GPSは複数の衛星までの距離と衛星位置を用いてユーザー位置を求める
- まず衛星までの距離と各衛星の位置を計算しなければならない
座標系: ECEF, WGS 84, ENU
- 地球に固定されてともに回転する直交座標系は ECEF(Earth Centered, Earth Fixed) と呼ばれる
- 地表に静止したユーザーの座標が時間とともに一定に保たれるため、ユーザー位置の表現に適している
- 最も一般的に使われるECEF座標系は、米国防総省が開発した WGS 84
- 原点は地球の質量中心
- z軸はCTP(Conventional Terrestrial Pole)を通る
- CTPは1900〜1905年の地球極位置の平均で、実際の極位置は約15m半径の円内を動く
- x軸はCTP赤道面と基準子午線であるMean Greenwich Meridianの交点を通る
- 衛星運動はニュートン力学に従って慣性座標系で扱うのが自然だが、GPS Interface Specificationは特定時点の衛星位置を ECEFフレーム で計算する手順を提供している
- 局所的な応用では、ユーザー位置を原点とする ENU(East-North-Up) 座標系が便利
- ECEF座標はユーザーの緯度・経度を用いた行列積でENU座標へ変換できる
- 衛星の方位角と仰角の計算にENU変換が使われる
高度の定義: 基準楕円体とジオイド
- 高度は「何を基準にするか」を先に定める必要がある
- 基準楕円体(reference ellipsoid) は地球を扁平な楕円体として抽象化したモデル
- 地球中心に位置し、回転軸はECEF z軸と一致する
- 地球を球として見るときによく使われる半径6371kmは、長半径と短半径の中間値
- 実際の地表上の点は基準楕円体の上にも下にもありうる
- ジオイド(geoid) は同じ重力ポテンシャル値を持つ点の集合であり、物理的意味を持つ面
- ジオイド基準の高度は正標高、または平均海面(MSL)上の高さと呼ばれる
- ジオイドは通常、基準楕円体上の高さ値の集合として指定される
- 緯度・経度・高度は楕円体座標として定義される
- 測地緯度は、点Pで楕円体表面に垂直な線と赤道面のなす角
- 地球中心と点Pを結ぶ線の角度は地心緯度であり、地球が完全な球なら測地緯度と一致する
- 楕円体座標から直交座標への変換は1段階で可能だが、ECEFから楕円体座標への変換には高速に収束する反復手続きが必要
衛星位置計算
- 理想的な衛星軌道は、6つの ケプラー軌道要素 で記述される楕円軌道
- 5つの要素は楕円の大きさ・形状と軌道面の向きを決める
- 6つ目の要素は特定のepochにおける衛星位置を決める
- 実際の衛星軌道は、地球組成の不均一性や太陽・月の重力の影響により完全な楕円ではない
- GPSはこれらの摂動を補正するため 16個の軌道パラメータ を放送している
- GPS Interface Specificationのtable 20-IVに、軌道補正を含む衛星位置計算手順がある
- ユーザー位置は受信時刻tで計算されるが、GPS信号は衛星からt-τ時刻に出発している
- 衛星位置は信号送信時刻t-τで計算する
- 信号が伝搬するτの間に地球が回転するため、衛星位置ベクトルを地球回転量だけ回転させて、受信時刻tのユーザーECEFフレームに合わせる
- これは衛星位置を単純に時刻tで計算するのとは同じではない
疑似距離と時計バイアス
- GPS受信機は衛星信号に含まれる送信タイムスタンプと受信機時刻を比較し、その差に光速を掛けて衛星までの距離を概算する
- この測定値が 疑似距離(pseudorange)
- 衛星時計と受信機時計が完全に同期し、信号が真空中を直線的に光速で伝搬するなら実距離と等しい
- 実際には時計オフセットや大気遅延のため、実距離とは異なる
- 衛星時計バイアスは位置誤差が数千メートルに達しうるため、必ず補正しなければならない
- GPS ephemerisメッセージの係数を用いた多項式と相対論項で計算する
- 多項式が補正の大部分を担い、相対論効果は衛星位置に応じて約1〜10m寄与する
- 受信機時計バイアスはユーザー位置とともに推定すべき未知数
- アルゴリズムでは時計バイアスを光速と掛けて距離単位で扱う
- 大気遅延は電離層成分と対流圏成分に分かれる
- 電離層遅延は一般に約25mの位置誤差を引き起こす
- 対流圏遅延は一般に約2mの位置誤差を引き起こす
- 本文の実験解析ではこれらの遅延を無視している
ユーザー位置と時計バイアス推定アルゴリズム
- 補正済み疑似距離測定値は、実際のユーザー-衛星距離、受信機時計バイアス、モデル化されていない誤差の和として表される
- ユーザー位置と時計バイアスは、測定疑似距離と予測疑似距離の差を最小化する値として求める
- 解法は 反復最小二乗 手続き
- ユーザー位置の初期値は
[0 0 0] - ユーザー時計バイアスの初期値は
0 - 各反復で現在の位置推定値に基づく衛星方向の単位ベクトルを積み上げてG行列を構成する
- 位置補正量と時計バイアス補正量を解き、変化量が閾値より小さくなるまで繰り返す
- ユーザー位置の初期値は
- 衛星がちょうど4個で非退化配置なら直接解を求められる
- 空が遮られていなければさらに多くの衛星が見え、一般には最小二乗解を使う
- 実装手順は次の流れに従う
- 生の疑似距離と衛星ephemerisを入力として受け取る
- 衛星ごとの時計バイアスを計算して疑似距離を補正する
- 可能なら電離層・対流圏補正を適用する
- 現在の受信機時計バイアスで疑似距離を補正する
- 疑似距離を光速で割って信号伝搬時間τを求める
- t-τ時刻の衛星位置を計算する
- τの間の地球回転を反映して衛星位置をユーザーECEFフレームに合わせる
- G行列と疑似距離差を作り、位置と時計バイアス補正量を計算する
Matlab実装上の詳細
- Matlabコードの大半は、右辺の既知の値から左辺の未知数を一度に評価する形になっている
- 一部の計算には解析的な閉形式解がなく、solver が必要
- 例は衛星位置計算中、平均近点角Mから離心近点角Eを求める段階
E - e*sin(E) == Mの関係は閉形式で解けないためvpasolveを使う
- 付録コードには次の機能が含まれる
- ユーザー位置と時計バイアスの計算
- 衛星位置の計算
- ユーザー位置と時計バイアスの最小二乗解計算
- 衛星時計バイアスの計算
- ECEF WGS84座標を楕円体座標へ変換
- ephemerisデータ形式の変換
生GPSデータ収集の構成
- 生のGPSデータを得るには、内部で位置だけ計算して出力する一般的なGPS機器ではなく、生の疑似距離や衛星ephemerisのような timing情報 を出力する受信機が必要
- u-blox NEO-M8T と 6T チップが要件に合う
- GPSユニット、アンテナ、シリアル出力ポートを含むハードウェア組み立て品をAmazonで約40ドルで購入できる
- 生GPS信号の受信と保存にはRTKLibの STRSVR ユーティリティを使用する
- RTKLibはGPS、Glonass、Galileo、BaiduなどGNSSの標準・精密測位を支援するオープンソースのプログラムパッケージ
- STRSVRはu-blox受信機のカスタム形式出力をRTCM標準形式へ変換する
- 必要な情報はRTCMメッセージ 1002 と 1019 に入っている
- 1002は生の疑似距離情報を含む
- 1019は衛星ephemeris情報を含む
- STRSVRはシリアルポート9600 Baudでデータを受け取り、RTCM 3形式ファイルとして保存するよう構成する
- データ収集はアパート建物の屋上で行われた
- GPS受信機は空が遮られない位置に置く
- u-blox u-centerソフトウェアで十分な衛星が見え、良好なposition fixが可能か確認する
- 約1時間の生GPSデータを収集する
RTCM処理とgoGPSの使用
- STRSVRは生GPSデータをバイナリ RTCM3 形式で保存する
- Matlabで処理するには、RTCM3データをデコードしてMatlabのデータ構造にする必要がある
- RTCMデコーダを自作する代わりに、goGPS Matlabライブラリの
load_stream関数を使用する- RTCM形式ファイルを読み、RTCMメッセージを抽出する
- 抽出したデータを
.matファイルとして保存し、位置計算アルゴリズムの入力に使う
rtcm_dataファイルも提供される- WordPressのセキュリティ制限のため
.matの代わりに.txt拡張子で提供される - ダウンロード後に
.matへ再度名前を変更する必要がある
- WordPressのセキュリティ制限のため
実験結果: 位置変動と時計ドリフト
- 受信機はデータ収集中は静止していたため、計算された位置の時間変化は位置計算アルゴリズムの実際の性能を示す
- ユーザー中心ENUフレームでの位置成分の標準偏差は次のとおり
- East: 14.00m
- North: 39.88m
- Up: 47.35m
- 位置変動はEastとNorth方向で約30m水準で、Up方向はより大きい
- 受信機時計バイアスは定数ではなく、時間とともに線形にドリフトする
- 時計バイアスはアルゴリズムでは距離単位で扱われる
- 結果プロットでは光速で割って時間単位へ変換される
- ドリフト量は 4.27e-7sec/sec
衛星方位角と仰角の計算
- 衛星の方位角と仰角はユーザー視点で定義されるため、ユーザー中心 ENUフレーム で計算する
- 計算手順は次のとおり
- ECEFフレームでユーザーから衛星へ向かう位置ベクトルを計算する
- ユーザー位置を楕円体座標である緯度・経度へ変換する
- その位置ベクトルをユーザー中心ENUフレームへ回転する
- ENU座標から方位角と仰角を計算する
- 例のepochで計算された8個の衛星の仰角はすべて正
- 方位角は正負の両方を取りうる
- ユーザーは地平線下の衛星を見ることができないため、仰角が正であるのは自然
- 複数のepochで同じ手順により衛星位置を計算すれば、GPS処理ソフトウェアが示す satellite track chart を作成できる
DOP: 位置推定品質の幾何学的要素
- DOP(Dilution of Precision) は位置推定がどれほど良いかを評価する指標
- 位置誤差には測定雑音だけでなく、ユーザー-衛星の 幾何配置 が影響する
- 疑似距離と衛星位置の測定がよりノイジーであるほど位置誤差は大きくなる
- 衛星が方位角・仰角で広く散らばっているほど幾何配置は有利で、DOPは低くなる
- 位置と時計バイアス誤差の共分散は、ユーザー距離誤差とG行列の関数として分解される
- G行列はユーザーから衛星へ向かう単位ベクトルで構成される
- ECEFフレームのG行列はDOP計算の便宜のためENUフレームへ回転する
- DOP成分はEast、North、Up方向で定義される
- HDOPはEastとNorth成分を結合した水平DOP
- VDOPはUp成分の鉛直DOP
- 実データではHDOPとVDOPは概ね 2.5未満
- この値は十分な水準とみなされる
- VDOPはHDOPより大きい
- 地表のユーザーは地平線下の衛星を観測できず、地平線上10度未満の衛星信号もノイズが大きすぎて通常使わないため、VDOPは高くなる
GPSインフラの規模
- GPSコンステレーションの構築には約 300億ドル がかかり、米国政府は維持に年間約 10億ドル を支出している
- GPSがなければ存在しえなかったUberの価値は700億ドル以上と言及されている
- GPSが可能にした多くの応用まで含めれば、GPSへの公共投資は非常に大きな経済的・技術的波及効果を生んだ事例と見なせる
1件のコメント
Hacker Newsのコメント
Androidはかなり前から搬送波位相（carrier phase）へのアクセスを提供しており、これを使うと、同じ近所くらいにあるデバイス同士の相対位置を、GNSSアンテナがデバイスのどこに隠れているかを気にしなければならないほどの精度で求められる
これだけではものすごく大したものではないが、各デバイスの加速度計やジャイロスコープまで組み合わせると良くなる
GNSSの擬似距離測定値の変化は、静止していなくても予測可能なので、リアルタイムで動作しながら性能低下も小さい
例えば、車輪のない模型飛行機をトラックの荷台に自動着陸させ、擦り傷や芝の滑走路への依存も避けられる
消費電力が非常に重要でなければ、かなり良いGNSS受信機も高価にする必要はないのに、なぜ100ドルでペアを普通に買えないのか分からない
ここに30ドルのGNSS受信機がある。不要なガラクタも少し付いてくるけど
自分でGPS受信機を作りたい人向けに、理論も多く説明している完全なオープンソースプロジェクトがある: http://www.aholme.co.uk/GPS/Main.htm
https://lea.hamradio.si/~s53mv/navsats/theory.html
GPSのコメントでよく挙がるが、それだけの理由がある記事: https://ciechanow.ski/gps/
別の解説で、おそらくよりインタラクティブ:
https://ciechanow.ski/gps/
もう一つ良いオープンソース実装がある:
https://m.youtube.com/watch?v=dVD1Yws__v0
ときどき、それもごく短時間だけ水面に上がる水生生物のGPSデータを集めている研究者を見たことがある
生データを記録して後処理すれば、消費電力と衛星信号にさらされる必要がある最小時間の両方を大幅に減らせ、露出時間は1秒未満まで下げられる
「下の図は、ユーザーとソースの幾何配置がユーザー位置の不確実性にどう影響するかを示している」という部分を見て、スマートフォンの地図アプリで位置の不確実性の形を円ではなく、こういう弧の交差形状に変える設定があればいいと思った
GPSは日常生活で相対論的効果を考慮しなければならない数少ない応用の一つだと聞いた。だとすると、生成されたデータにはすでにこうした相対論的効果が取り除かれているということで合っている？
商用GPSデバイスの出力なら、その通り。出力を生成するための取得後処理で、さまざまな誤差要因が補正されている
この記事は複数の衛星からストリーミングされる生のGPSデータを扱っており、そこから出力値を作るには処理が必要で、精度を高めるために地上局や海上補正のような追加入力が入る場合も多い
複数のGPS機器ベンダーは概ね似たようなことをしているが、細部が肝心
他のコメントがリンクしている https://ciechanow.ski/gps/ も読みやすい
観測者より速く動く物体では時間もより遅く流れ、衛星はかなり速く動いている
GPSは観測者と衛星の間で時間が同期していなければならないため、特殊相対性理論と一般相対性理論の効果を反映するよう時間の流れが補正されている
そのため相対論は重要だが、自分の位置を解くのに相対論を深く知る必要はない
ただし、さまざまな形の長基線／ネットワークRTKでは、より精緻なモデリングが必要になることがある
次のステップはPPPやRTK。GNSSは可能性の沼がとても楽しい
地球平面論者向けの練習問題: 球状の地球を周回する衛星なしに、スマートフォンのGPS地図がどう動作するのか説明せよ。解答過程を示せ
まず、衛星なしではスマートフォンがそれを実現できないという、かなり複雑な論証をしなければならない
それより少し簡単で、地球平面論者がかわしにくい例としては、ISSは肉眼でもほとんど見え、裏庭の望遠鏡なら確実に見えるという点がある。Starlink衛星も同様
地球平面論は理性によって到達した立場ではなく、ほとんど常に混乱から生じるか、揺るがない中核的信念の必然的な帰結として生じる。たいていは聖書を極端に文字通り読むこと、あるいは「公式なものはすべて嘘」という偏執的な妄想に由来する