s-taka.org

高橋 賢

Lang:

QZS L6 ToolによるMADOCA-PPP電離層遅延情報の表示

category: gnss
tags: madoca-ppp qzss rtklib

はじめに

内閣府が運用する準天頂衛星（QZS: quasi-zenith satellite）みちびきは、GPSのような測位信号を放送するだけでなく、測位精度を高める信号も放送しています。

この信号のひとつに、L6周波数帯にて放送されるMADOCA-PPP（マドカ・ピーピーピー、multi-GNSS advanced orbit and clock augmentation - precise point positioning）があります。将来のみちびき衛星からは、MADOCA-PPPに新しく電離層遅延情報が放送される予定です。

また、内閣府がオープンソースとしてGitHubに公開しているMADOCALIB（マドカリブ、MADOCA-PPP test library）では、バージョン1.2からこの電離層遅延情報を扱えます。この電離層遅延情報により、短時間にて高精度測位ができました。

今回、この電離層遅延情報の内容を観測するために、オープンソースQZS L6 Tool（キューズイエス・エルシックス・ツール）をアップデートしました。

MADOCA-PPPの電離層遅延情報

みちびき衛星は、L6D信号にて日本国内向けのCLAS（シーラス、centimeter level augmentation service）と、L6E信号にて世界で利用できるMADOCA-PPPを放送しています。今後打ち上げられる新しいみちびき衛星から、L6D信号でMADOCA-PPPの電離層遅延情報が伝達される見込みです。この信号フォーマットは、仕様書IS-QZSS-MDC-002に記述されています。

しかしながら仕様書を読んでも、電離層遅延情報がどのように伝送されるか、わかりませんでした。そこで、MADOCA-PPPアーカイブデータを用いて、伝送方法を確認してみます。アーカイブデータには、将来のみちびき衛星から放送されるであろうデータも先行してアーカイブされています。ここでの対象PRN（pseudo random noise）番号は200番と201番の2つです。電離層遅延情報の提供開始は2024年6月28日00:00 UTCです。

CLASやMADOCA-PPPは、2000ビットをデータパート（DP）として、1秒間に1DPを伝送します。DPにはヘッダやリードソロモン誤り訂正符号が含まれ、これらを除いた正味のデータ長は1695ビットです。そして5DPで1サブフレーム（SF）を、6SFで1フレームを構成します。すなわち、これらの情報は30秒周期で伝送されます。伝送されるデータの種別は、サブタイプ（ST）番号にて識別されます。このSTデータは可変長で、SF内において前に詰めて伝送されます。一つのDPにて複数STデータが含まれることや、複数DPにまたがって一つのSTデータ伝送されることもありますが、SF間をまたがってSTデータが伝送されることはありません。SF終端部分の余りビット列はゼロパディングされます。SF先頭となるDPのヘッダには、1ビット長のsubframe indicatorがオンになります。また、DP先頭にある12ビット長のmessage numberには、4073の固定値が用いられます。

一方、新しいMADOCA-PPP電離層遅延情報には、DP先頭にある12ビット長のmessage type（MT, 上述のmessage numberと同義）として、1または2が用いられます。MT1は座標定義のためのメッセージであり、1メッセージに複数の長方形または円形の平面がエリア番号とともに定義されます。一方、MT2では、MT1で定義された平面内での衛星ごとの擬似距離修正量が伝達されます。これらも可変長データです。最初に、subframe indicatorオンのDPにて、MT1メッセージが一つ伝達され、定義された平面エリア番号数だけの複数MT2が続きます。MT1は、8ビット長の地域ID（region ID）ごとに伝送されますが、地域IDと座標との関係はわかりません。複数DPにまたがって情報伝送されることもあり、subframe indicatorオフのDPは、subframe indicatorオンのDPに連結させなければなりません。現時点では、連結される最大DP数は定義されていません。余りビット列については、CLASなどと同様に、ゼロパディングされます。一方、CLASなどと異なり、送るべき情報がないときには、subframe indicatorオンのオールゼロ内容のDPが伝送されます。CLAS、MADOCA-PPP、MADOCA-PPP電離層遅延の伝送方法を表にまとめます。

CLASMADOCA-PPPMADOCA-PPP電離層遅延情報
PRN番号範囲193-199203-211200, 201
データ識別ST番号（MT番号は4073固定）同左MT番号
データ周期30秒間同左不定
SF先頭内容ST1またはST3同左MT1

QZS L6 ToolによるMADOCA-PPP電離層遅延情報の表示

電離層遅延情報の伝送方法がわかりましたので、これをもとに電離層遅延情報を復号するコードを作成します。打ち上げ予定のみちびき5号機、6号機、7号機から放送されるであろうL6データも先行配信するMADOCA-PPP補強情報のインターネット・ストリーム配信を利用して、リアルタイム表示してみました。アーカイブデータは誰でも利用できますが、このストリーム配信の利用には申請が必要です。

次の映像の左がPRN200番のもの、右がPRN201番のものです。MADOCA-PPPでは、1秒間に1DPが伝送されますので、この映像も1秒ごとに変化します。

ここでは、トレースオプション-t 1にて、補強情報の内容を表示しています。具体的には、str2str -in ntrip://[user id]:[password]@[NTRIP caster address]/[mount point] 2> /dev/null | qzsl6read.py -t 1などとしています。ヌル以外の情報が伝達されたときには多くの文字列が表示されます。この映像から、PRN200とPRN201とで異なる情報が伝送されていることがわかります。また、ヌルDPも多く観測されています。

メッセージタイプ1による平面座標定義

ここで、一例として、2024年8月16日 00:00:00 UTCから1時間のPRN 200のデータをみちびきアーカイブサイトからダウンロードし、PRN200のMT1情報を観測してみます。上述の文字列からテキストエディタやページャーにて先頭にMT1が含まれる文字列を抜き出します。

MT1 Epoch=00:00:30+5 UI=30s(5) MMI=0 IODSSR=2 Region=1  4191bit  NumAreas=3
 # shape lat[deg] lon[deg] lats lons / radius[km]
 5 RECT     -34.0    132.5  3.0  6.5
 6 RECT     -31.5    119.5  4.0  6.5
 8 RECT     -20.5    119.5  7.0  6.5
MT1 Epoch=00:00:33+5 UI=30s(5) MMI=0 IODSSR=2 Region=3  2504bit  NumAreas=2
 # shape lat[deg] lon[deg] lats lons / radius[km]
 1 RECT      17.4    121.4  1.3  1.1
 2 RECT      14.7    120.8  1.4  1.0
MT1 Epoch=00:00:35+5 UI=30s(5) MMI=0 IODSSR=0 Region=4* 77bit  NumAreas=1
 # shape lat[deg] lon[deg] lats lons / radius[km]
 1 RECT      -6.8    107.0  0.9  1.2

ここでは、地域ID 1に3長方形平面が定義されています。最初のデータは、エリア番号5の長方形平面で、南緯34.0度と東経132.2度（オーストラリア・アデレード西部の海中）を基準として、緯度±3度、経度±6.5度の幅（span）が定義されています（一般にspanとは先端から先端までの距離を指しますが、ここでは中央から先端までの角度差を表します）。同様に、地域ID 3に2平面、地域ID 4に1平面がそれぞれ定義されています。ここでは、地域ID 4にアラート（アスタリスク*で表されています）が出ていました。そして、この3地域IDを含むMT1が、繰り返し伝送されていました。

一方、PRN 201 MT1の内容は次の通りです。

MT1 Epoch=00:00:30+5 UI=30s(5) MMI=0 IODSSR=2 Region=2  2614bit  NumAreas=2
 # shape lat[deg] lon[deg] lats lons / radius[km]
 4 RECT     -37.5    146.5  6.5  7.5
 7 RECT     -26.0    132.5  5.0  6.5
MT1 Epoch=00:00:32+5 UI=30s(5) MMI=0 IODSSR=0 Region=5  9936bit  NumAreas=8
 # shape lat[deg] lon[deg] lats lons / radius[km]
 1 RECT      43.7    142.5  2.5  4.3
 2 RECT      39.2    140.5  2.0  3.0
 3 RECT      34.2    140.0  3.0  2.5
 4 RECT      35.2    135.0  3.3  2.5
 5 RECT      32.8    130.3  2.8  2.2
 6 RECT      27.5    129.0  2.5  3.0
 7 RECT      24.5    124.2  1.0  1.8
 8 CIRCLE    26.9    142.2  100

ここでは、地域ID2と5が繰り返して定義され、それぞれに2平面と8平面が定義されていました。地域ID 5、エリア番号8には、小笠原の父島と母島との中間にある海上を中心とした半径100キロメートルの円形平面が定義されていました。

同様に、MADOCALIB ver.1.2に同梱されているPRN200のデータ2024162all.200.l6を読み出した結果は次のとおりです。これは、MADOCALIBのsample_data/にあります。サンプルデータは、北緯36.11度、東経140.09度近辺にある国土地理院のデータを処理するもので、2024年の1月1日からの日数が162日目（6月10日）のものです。

MT1 Epoch=00:00:30+1 UI=30s(5) MMI=0 IODSSR=2 Region=1  4551bit  NumAreas=3
 # shape lat[deg] lon[deg] lats lons / radius[km]
 5 RECT     -34.0    132.5  3.0  6.5
 6 RECT     -31.5    119.5  4.0  6.5
 8 RECT     -20.5    119.5  7.0  6.5
MT1 Epoch=00:00:33+1 UI=30s(5) MMI=0 IODSSR=2 Region=3  2154bit  NumAreas=2
 # shape lat[deg] lon[deg] lats lons / radius[km]
 1 RECT      17.4    121.4  1.3  1.1
 2 RECT      14.7    120.8  1.4  1.0
MT1 Epoch=00:00:35+1 UI=30s(5) MMI=0 IODSSR=0 Region=4  1077bit  NumAreas=1
 # shape lat[deg] lon[deg] lats lons / radius[km]
 1 RECT      -6.8    107.0  0.9  1.2

このPRN200 MT1の内容は、先の8月16日 00:00:00 UTCから1時間のPRN200 MT1の内容と同一でした。この平面座標定義列に上述の座標と一致するものはありませんでした。MADOCALILB ver.1.2で計算したPRN200のものでは、電離層遅延補強が効いていなかったようです。次に、PRN201のデータ2024162all.201.l6も表示してみます。

MT1 Epoch=00:00:30+1 UI=30s(5) MMI=0 IODSSR=2 Region=2  2794bit  NumAreas=2
 # shape lat[deg] lon[deg] lats lons / radius[km]
 4 RECT     -37.5    146.5  6.5  7.5
 7 RECT     -26.0    132.5  5.0  6.5
MT1 Epoch=00:00:32+1 UI=30s(5) MMI=0 IODSSR=0 Region=5  8536bit  NumAreas=8
 # shape lat[deg] lon[deg] lats lons / radius[km]
 1 RECT      43.7    142.5  2.5  4.3
 2 RECT      39.2    140.5  2.0  3.0
 3 RECT      34.2    140.0  3.0  2.5
 4 RECT      35.2    135.0  3.3  2.5
 5 RECT      32.8    130.3  2.8  2.2
 6 RECT      27.5    129.0  2.5  3.0
 7 RECT      24.5    124.2  1.0  1.8
 8 CIRCLE    26.9    142.2  100

このPRN201 MT1の内容もまた、8月16日のPRN201 MT1の内容と同一でした。上述の国土地理院の座標に一致している平面座標定義は、地域ID 5、エリア番号3にありました。北緯34.44度、東経132.41度にある広島市も、PRN201の地域ID 5、エリア番号5に属しています。

上述の私の平面座標表記はわかりにくいので、今後、QGISなどを用いてわかりやすい図面で表したいと思います。

現在のMADOCA-PPP電離層遅延情報に含まれる地域数は限られていますが、回線容量には多くの余裕があります。電離層遅延情報の今後にも期待しております。

1月1日からの経過日数計算

ところで、このような解析には、1月1日からの日数やGPS週番号をカレンダー形式で表示できるソフトウェアがあると便利です。私は、以前に国土地理院のウェブサイトに掲載されていたgpscalを利用しています。1996年7月20日のHatanaka Y.さんのクレジットのある、64行からなるc-shellスクリプトです。シェルスクリプトで日数を計算した上で、内部からcalコマンドを呼び出し、awkコマンドで整形しています。素晴らしいソフトウェアなので、再び公開していただけると嬉しいです。

$ gpscal.csh 06 2024
               June 2024
Week   Sun Mon Tue Wed Thu Fri Sat
2316                             1
                               153
2317     2   3   4   5   6   7   8
       154 155 156 157 158 159 160
2318     9  10  11  12  13  14  15
       161 162 163 164 165 166 167
2319    16  17  18  19  20  21  22
       168 169 170 171 172 173 174
2320    23  24  25  26  27  28  29
       175 176 177 178 179 180 181
2321    30
       182
$

セマンティック・バージョニング

QZS L6 Toolにセマンティック・バージョニングを採用することにしました。これは、major.minor.patchの3個の数字部分からなるバージョン番号付与方法です。出力に互換性を失うときにはmajor番号を、後方互換性を保ったままで機能追加するときはminor番号を、バグ修正するときはpatch番号を更新します。表示桁数の調整などは、patch番号更新にするつもりです。現在を元期にして、QZS L6 Toolのバージョン番号を0.1.0としました（従来のカレンダ・バージョニングでは20240816updです）。

まとめ

みちびきMADOCA-PPPの電離層遅延補強情報を復号するコードを作成しました。

電離層遅延情報は、異なる衛星でそれぞれPRN200とPRN201とを用いて伝送される予定です。先行して公表されているアーカイブデータやストリーミングを観測したところ、PRN200とPRN201とで異なる内容が伝送されていました。現在のところ、それぞれのPRNで対象地域が固定されているようですので、おおまかな座標から受信すべきPRNを決められるようです。例えば、北緯36.11度、東経140.09度近辺にある国土地理院では、みちびき7号機から放送されるであろうPRN201を利用することになります。

今後、解析を進め、効果的なMADOCA-PPPの利用方法を考えてみたいと思います。

関連記事