Galileo HAS Reference User Algorithm in Matlab
はじめに
アメリカのGPS、ヨーロッパ連合のGalileo、日本の準天頂衛星みちびきなど、測位衛星システム（GNSS: Global Navigation Satellite System）が放送する電波を受信することにより、私たちは自らの座標を推定できるようになりました。この電波受信から得た衛星座標、時刻、衛星から受信機に至る擬似距離に対し、補正により測位精度を高める方法は補強（オーギュメンテーション）と呼ばれます。
Galileoでは、HAS（High Accuracy Service）という補強サービスを実施しています。この補強情報は、Galileo衛星のE6-B信号にて伝達されるほか、インターネットでも配信されています。この補強情報は、誰でも無料で利用できます。
このHASについて、性能評価評価ソフトウェアとサンプルデータが、欧州のGSC（European GNSS Service Center）より一般に提供されましたので、試してみました。
ソフトウェアとサンプルデータの入手
このソフトウェアの利用には、Matlabと、そのStatistics and Machine Learning Toolboxが必要です。Windows環境下では、さらにParallel Computing Toolboxの導入により、計算が並列化でき、処理時間を短縮できます。
GSCのメーリングリストにて、このソフトウェアとデータの提供開始を知り、2026年6月20日に申請しました。
この申請には、あらかじめGSCアカウントの登録が必要です。登録が済んだら、ログインし、この「Request access」ボタンをクリックします。
遷移した利用許諾画面を読み、最後に、氏名と電子メールアドレスの入力と、国名と利用許諾の選択をした上で「Submit」ボタンを押します。
次の画面では、ZIPファイルへのダウンロードリンクが表示されますので、ダウンロードします。容量は約3.5ギガバイトあり、パスワードロックされています。私の環境では、ダウンロードし終えるまでに半日かかりました。ZIPファイルのパスワードは、約1週間後に電子メールで伝えられました。
パスワードを入手し、ZIPファイルを解凍すると、そこには複数のZIPファイルと、Windows用解凍スクリプトuncompress.cmd、そして、ソフトウェアの利用方法が書かれたファイルがありました。
ここにあるZIPファイル全てを、ディレクトリー構造を保ったままで解凍する必要があります。私は、最初、Mac上の解凍ソフトウェアKekaで解凍を試みましたがディレクトリー構造を保つことが難しく、結局、Windows PC上のコマンドプロンプトで、このuncompress.cmdを実行して解凍しました。解凍後の容量は、全部で7ギガバイト程度あります。ソースコードは、MatlabのMファイル（テキスト形式）ですので、Windows上でもMac上でも同様に動作します。
ディレクトリー構造
トップディレクトリーには、以下の2つのMファイルと、4つのサブディレクトリーconf、data、doc、srcがありました。
- Example_Run_HASUA_PPPSD_stat_year_doy_all.m （計算プログラム）
- Example_Plot_HASUA_PPPSD_stat_year_doy_any.m （結果表示プログラム）
サンプルデータは、dataディレクトリーの中にあります。ここには、地点別のRINEX形式の観測データを収めたディレクトリー（地点名BADG、MET3、SEYG、THU2、VILL）のほかに、RINEX形式航法データや補強データやアンテナ情報ANTEXを収めたcommonディレクトリがあります。
それぞれのディレクトリーは、年番号（ここでは2023年のみ）と、DOY（day of yearのことで1月1日を1日目とした経過日数、ここでは356から365）に分類されています。計算結果は、これらのディレクトリーの中のresults_で始まるディレクトリに収められます。
計算は、Galileo E1-E5a, GPS L1-L2、E1-E5a, L1-L2、E1-E6, L1-L2の3通りの条件に加え、インターネット配信IDD(Internet Data Distribution)とE6-B信号による配信SIS(Signal-in-Space)との2通りの、合計6通りについて行われます。
計算結果として、計算実行日時のついたファイル名にて、MAT形式の結果ファイル、拡張子solのテキスト形式の結果ファイル、および拡張子logのログファイル、の3種類が出力されます。この結果ディレクトリには、あらかじめ計算された結果ファイルも（例えばBADG00RUS_R_20233560000_01D_30S_MO_PPPSD_160120_SSRC_SIS_2026_04_28_12_59_12.mat）収められていました。
data
├── BADG
│ └── 2023
│ ├── 356
│ │ ├── results_E1_E5a_L1_L2_IDD
│ │ ├── results_E1_E5a_L1_L2_SIS
│ │ ├── results_E1_E5b_L1_L2_IDD
│ │ ├── results_E1_E5b_L1_L2_SIS
│ │ ├── results_E1_E6_L1_L2_IDD
│ │ └── results_E1_E6_L1_L2_SIS
│ ├── 357
...
├── common
│ └── 2023
│ ├── 356
│ ├── 357
│ ├── 358
│ ├── 359
│ ├── 360
│ ├── 361
│ ├── 362
│ ├── 363
│ ├── 364
│ └── 365
├── MET3
...
ソースコードはsrcディレクトリーに収められています。MatlabのMファイルは、1関数を1ファイルにて記述します。関数の数を数え上げてみると383にもなり、大きなプログラムであることがわかります。
du -a src/ | grep '\.m' | wc -l
383
入力データ
commonディレクトリーにある合計7日分のRINEX形式の航法データのうち、最初のもの（DOY=356、2023-12-22）の先頭部分は次のとおりです。この航法データには、日本の準天頂衛星みちびきやインドのNavICのものも含まれていました。
3.05 NAVIGATION DATA MIXED RINEX VERSION / TYPE
BCEmerge congo 20231222 004605 GMT PGM / RUN BY / DATE
gfzrnx-2.1.9 FILE MERGE 20241112 113001 UTC COMMENT
BDSA 4.1910e-08 3.7253e-08 -1.0133e-06 1.6093e-06 IONOSPHERIC CORR
BDSB 1.1469e+05 1.8022e+05 -1.8350e+06 1.9661e+06 IONOSPHERIC CORR
GAL 1.4300e+02 -8.9844e-01 5.5237e-03 0 IONOSPHERIC CORR
GPSA 2.6077e-08 7.4506e-09 -1.1921e-07 1.1921e-07 IONOSPHERIC CORR
GPSB 1.4950e+05 -2.1299e+05 0.0000e+00 3.2768e+05 IONOSPHERIC CORR
IRNA 8.1025e-08 2.9802e-07 -2.8610e-06 -7.5102e-06 IONOSPHERIC CORR
IRNB 1.2493e+05 7.3728e+05 -2.0972e+06 -8.3231e+06 IONOSPHERIC CORR
QZSA 6.6124e-08 -7.0035e-07 1.9670e-06 0.0000e+00 IONOSPHERIC CORR
QZSB 7.9872e+04 1.1960e+06 -8.3886e+06 -8.3886e+06 IONOSPHERIC CORR
GAGP 9.6042640507e-10 8.881784197e-16 432000 2293 TIME SYSTEM CORR
GAUT -1.8626451492e-09 8.881784197e-16 345600 2293 TIME SYSTEM CORR
GLGP -3.9115548134e-08 0.000000000e+00 345600 2293 TIME SYSTEM CORR
GLUT 9.3132257462e-10 0.000000000e+00 345600 2293 TIME SYSTEM CORR
GPUT 0.0000000000e+00 4.440892099e-15 61440 2294 TIME SYSTEM CORR
IRGL 5.5530108511e-08-4.440892099e-14 430800 2293 TIME SYSTEM CORR
IRGP 4.5984052122e-09-1.332267630e-15 430800 2293 TIME SYSTEM CORR
QZUT 4.6566128731e-09 0.000000000e+00 8192 2294 TIME SYSTEM CORR
18 18 1929 7 LEAP SECONDS
また、一例として、地点BADGの同一日の観測データの先頭部分は次のとおりです。ここから、測位は30秒ごとに行われ、24時間分のデータであることがわかります。HASの補強対象はGalileoとGPSですが、観測データにはロシアのGLONASSのものも含まれていました。
3.04 OBSERVATION DATA M RINEX VERSION / TYPE
JPS2RIN v.2.0.178 JAVAD GNSS 20231223 000638 UTC PGM / RUN BY / DATE
AUTOMATIC IAA OBSERVER / AGENCY
BADG MARKER NAME
12338M002 MARKER NUMBER
02682 JAVAD TRE_3 DELTA 3.7.10 Oct,22,2020 REC # / TYPE / VERS
-838282.9631 3865774.0774 4987620.6604 APPROX POSITION XYZ
00328 JAVRINGANT_DM JVDM ANT # / TYPE
0.0280 0.0000 0.0000 ANTENNA: DELTA H/E/N
G 20 C1C L1C D1C S1C C1W L1W D1W S1W C2X L2X D2X S2X C2W SYS / # / OBS TYPES
L2W D2W S2W C5X L5X D5X S5X SYS / # / OBS TYPES
R 20 C1C L1C D1C S1C C1P L1P D1P S1P C2C L2C D2C S2C C2P SYS / # / OBS TYPES
L2P D2P S2P C3X L3X D3X S3X SYS / # / OBS TYPES
E 20 C1X L1X D1X S1X C8X L8X D8X S8X C6X L6X D6X S6X C7X SYS / # / OBS TYPES
L7X D7X S7X C5X L5X D5X S5X SYS / # / OBS TYPES
26 R01 1 R02 -4 R03 5 R04 6 R05 1 R06 -4 R07 5 R08 6 GLONASS SLOT / FRQ #
R09 -2 R10 -7 R11 0 R12 -1 R13 -2 R14 -7 R15 0 R16 -1 GLONASS SLOT / FRQ #
R17 4 R18 -3 R19 3 R20 2 R21 4 R22 -3 R23 3 R24 2 GLONASS SLOT / FRQ #
R25 -5 R26 -6 GLONASS SLOT / FRQ #
30.000 INTERVAL
2023 12 22 0 0 0.0000000 GPS TIME OF FIRST OBS
2023 12 22 23 59 30.0000000 GPS TIME OF LAST OBS
18 LEAP SECONDS
高精度測位の正解である、それぞれの観測地点座標は、結果表示プログラムExample_Plot_HASUA_PPPSD_stat_year_doy_any.mに、ECEF（Earth-Centered, Earth-Fix）形式で書かれていました。これらをQZS L6 Toolのecef2llh.pyにて緯度・経度・楕円体高に変換して、以下にまとめます。
$ ecef2llh.py -.838282176007901E+06 0.386577732440940E+07 0.498762455124815E+07
51.7697041 102.2349908 811.421
| station name | latitute [deg] | longitude [deg] | ellipsoidal height [m] | Google map link |
|---|---|---|---|---|
| BADG | 51.7697041 | 102.2349908 | 811.421 | Badary RTF-32, Russia |
| MET3 | 60.2174569 | 24.3945036 | 79.232 | Metsähovi Geodetic Research Station, Finland |
| SEYG | -4.6787305 | 55.5306332 | -37.619 | Seycelles International Airport? |
| TASH | 41.3280498 | 69.2955723 | 439.713 | Teleskop, Uzbekistan |
| THU2 | 76.5370484 | -68.8250532 | 36.238 | Dundas?, Greenland |
| VILL | 40.4435961 | -3.9519751 | 647.341 | European Space Agency, Spain |
HAS計算と結果プロット
計算は、トップディレクトリーにあるExample_Run_HASUA_PPPSD_stat_year_doy_all.mを使います。観測データは、6地点、10日分、6通りの信号と配信経路からなるので、合計360種あります。このプログラムを実行すると、それぞれの観測データで、30秒ごとの1日分（2880エポック）の測位計算が行われることになります。
私は、Apple Silicon M1で16GBメモリー内蔵のMacデスクトップPCや、Core i3-8100で32GBメモリー内蔵のWindowsデスクトップPCなどで試しましたが、いずれも一つの観測データ処理につき、360秒程度の処理時間を要しました。私のWindows上MatlabでのParallel Computing Toolboxでは2並列にしかならず、また、Mac上にこのToolboxがあると、memory関数利用不可を示すエラーが発生します。Parallel Computing ToolboxがインストールされたMac PCでは、それが利用されないよう、ソースコードにhas_parallel=falseを追加する必要がありました。
実は、私のPCでは、すべての計算を完了することはできませんでした。計算させたまま帰宅し、翌日に出勤したときにはリブートしていたり、または、画面表示が壊れていて、計算が止まっていました。
計算結果のプロットには、Example_Plot_HASUA_PPPSD_stat_year_doy_any.m を用います。実行すると、地点名とDOYが問われ、結果ディレクトリーが表示されます。その中のMATファイルをクリックすると、30秒ほどで、以下のような12枚のグラフが次々と表示されます。これは、地点BADGの2023年12月22日のE6-B配信（SIS）の例です。
多くの専門的な解析結果が表示されました。今後、じっくりとこれらを読み、自らのデータでも試してみたいと思います。
まとめ
Galileo HASの性能評価ソフトウェアを試してみました。このソフトウェアはMatlabで書かれていて、ソースコードを読むこともできます。添付されているデータは、10日分で6地点からなりました。今後、自らのデータでも試してみたいと思います。
関連記事
- Galileo HAS（high accuracy service）その3 27th October 2023
- Galileo HAS（high accuracy service）その2 20th February 2023
- Galileo HAS（high accuracy service）その1 7th February 2023