s-taka.org

高橋 賢

Lang:

QZS L6 Toolの出力形式変更

category: gnss
tags: b2b clas mosaic qzsl6tool qzss

はじめに

測位衛星が放送するメッセージを表示するツールQZS L6 Toolを公開しています。

この名称は、準天頂衛星みちびき（QZS: quasi-zenith satellite）のL6周波数帯から名付けました。みちびきがL6周波数帯で放送するCLAS（シーラス、センチメーター級測位補強サービス、centimeter-level augmentation service）の高精度測位の仕組みを知りたかったからです。その後、欧州のGalileoのHAS（ハス、高精度サービス、high accuracy service）や、中国のBeiDouのPPP-B2b（precise point positioning - B2b-signal）のメッセージもQZS L6 Toolで読めるようになりました。

しかしながら、その表示形式が「項目=値」であったために、再利用や、システム間のメッセージの比較が困難でした。そこで、QZS L6 Tool出力の機械可読（machine-readable）と人間可読（human-readable）の両立を目指して、出力形式を変更しました。

QZS L6 Toolのコードを見直す過程で、CLASの仕様書を読み直し、バグも修正しました。仕様書記述を私が誤って解釈していた部分があることに気づく一方、現在でもわからないところをメモしました（おわかりの方がいらっしゃいましたら、教えてくださると嬉しいです！）。

ダウンロード

QZS L6 Toolは、ZIP形式でもダウンロードできます。しかしながら、gitを使ってダウンロードするのが、アップグレードの点でお勧めです。gitを用いてダウンロードするには、

git clone https://github.com/yoronneko/qzsl6tool

を実行します。すでにgitを用いてダウンロードされているのでしたら、そのディレクトリ内で、git pullを実行すればアップグレードされます。皆さんが独自の改良をされていても、gitがなるべく変更を追跡するように更新してくれ、また特定コードの変更箇所もわかりやすく表示できるので、便利です。

表示形式変更

みちびきCLASを例に、下のスクリーンショットでこの表示形式変更を比較します（左が変更前、右が変更後）。例えば、ST2（サブタイプ2、衛星軌道の修正量）の表示では、1行あたりの文字数が減り、読みやすくなったと思います。

QZS L6 Tool display style modification

新形式では、各データ列の最初に見出し（例えばST2 SAT IODE radial[m] along[m] cross[m]）が表示されます。衛星時刻補正では、clockとc0とが混在していましたが、両者とも同じものを指すので、c0に統一しました（衛星時刻補正にはc1やc2も用いられることがあります）。また、CLASでは、ネットワークIDだけでなくその地域名も、またグリッドではなく緯度・経度座標を表示するようにしました。

たったこれだけですが、資料をまとめるのには、便利です。例として、中国の測位衛星BeiDouが放送する高精度補強情報PPP-B2bのデータを処理してみます。QZS L6 Toolのsample/20230819-081730hasbds.septには、私がSeptentrio mosaic-go X5受信機で取得したデータがあります。これをtest/do_test.shにて処理した結果がtest/expected/20230819-081730hasbds.b2b.txtです。この内容の一部は次のとおりです。

C60 MT4  CLOCK 08:17:34 IODSSR=1 IODP=2
SAT IODCorr   c0[m]
C19       4   0.483
C20       4  -0.035
C23       6  -0.050
C60 MT4  CLOCK 08:17:34 IODSSR=1 IODP=2
SAT IODCorr   c0[m]
C24       6  -0.178
C26       4   0.000
C35       2   1.682
C37       3  -0.920
C39       2   0.339
C42       1   0.578
C45       0   0.443

これは、BeiDouの静止衛星C60が放送する時刻補正情報（MT4）であり、C19衛星やC20衛星などに対するものです。この1列目が衛星番号、2列目がデータ管理番号、3列目がメートル単位の衛星時刻補正量（時刻補正量に光速をかけて、距離単位にて表したもの）です。

この部分をコピーし、グラフ表示ソフトウェアgnuplotにて表示します。次のような内容のスクリプトファイル（例えばppp-b2b.gp）を作成し、先のデータを$data << EODとEODの間にペーストして（EODは別の文字列でもかまいません）、そのペースト部分について少々のフォーマット変更をします。

ppp-b2b.gp

# -*- coding: utf-8 -*-
# gnuplot script
# clock correction distance, source:
# directory of test/expect/20230819-081730hasbds.b2b.txt of
# QZS L6 Tool, https://github.com/yoronneko/qzsl6tool.

$data << EOD
# SAT IODCorr   c0[m]
C19       4   0.483
C20       4  -0.035
C23       6  -0.050
C24       6  -0.178
C26       4   0.000
C35       2   1.682
C37       3  -0.920
C39       2   0.339
C42       1   0.578
C45       0   0.443
EOD

set termoption enhanced
set termoption font "Times New Roman, 18"
set termoption linewidth 3
set border lw 0.5
set yrange[-2:2]; set ytics 1
set format y "%3.1f"
set xlabel "BDS satellite PRN"
set ylabel "clock correction distance [m]"
set key off
plot $data u (column(0)):3:xtic(1) ti col with boxes

# EOF

このフォーマット変更とは、見出し行の先頭のコメントアウトと、複数行にまたがるデータ間にある見出し行などの削除です。gnuplotでは、スクリプトファイル中の#で始まる行は、コメントアウトとみなされ、無視されます。

gnuplotを実行し、このスクリプトファイルをロードして、ボックスプロットします。load "ppp-b2b.gpにてグラフが表示され、set output "ppp-b2b.png" set term png repにてグラフがPNGファイルppp-b2b.pngに保存され、quitでgnuplotが終了します。

QZS L6 Tool display style modification

$ gnuplot -d

	G N U P L O T
	Version 6.0 patchlevel 0    last modified 2023-12-09

	Copyright (C) 1986-1993, 1998, 2004, 2007-2023
	Thomas Williams, Colin Kelley and many others

	gnuplot home:     http://www.gnuplot.info
	faq, bugs, etc:   type "help FAQ"
	immediate help:   type "help"  (plot window: hit 'h')

	Terminal type is now qt
gnuplot> load "ppp-b2b.gp"
gnuplot> set output "ppp-b2b.png"
gnuplot> set term png

Terminal type is now 'png'
Options are 'truecolor nocrop enhanced butt size 640,480 font "arial,12.0" '
gnuplot> rep
gnuplot> quit

もちろん、データをファイルに保存してプロットすることも、プロットファイル形式としてPDFを指定することも可能です。

グラフ作成後の、グラフの表示形式、軸名、縦横比などの変更、また、データ内容の差し替えや多数データのグラフ作成は、研究者や学生にとってよくある日常の出来事です。先生や上司に

  • すべてのグラフについて、縦軸と横軸の表示範囲を揃えなさい。
  • 数値間隔は1か2か5にしなさい。
  • すべてのグラフの作り直し！

と言われても、怒ってはいけません。グラフを再作成するには、gnuplotならば、スクリプトファイルを変更して、そのファイルをロードするだけです。多数のデータファイルを一括処理するバッチファイルを生成系AIに作ってもらっても良いでしょう。グラフ作成については、エクセルよりも、gnuplotの方が便利です。LaTeXにて報告書作成することにすれば、グラフ差し替えの完全自動化も可能です（私はそうしています！）。

gnuplotは、

  • Macのパッケージマネージャhomebrewならばbrew install gnuplotにて、
  • Windows標準搭載のパッケージマネージャwingetならばwinget install gnuplot.gnuplotにて、
  • Linuxディストリビューション標準搭載のパッケージマネージャaptならばapt install gnuplotにて、

簡単にインストールできます。

みちびきCLAS仕様書でわからなかったこと

CLASの仕様書(IS-QZSS-L6-005)、説明会資料、サンプルプログラム（CLASLIB）は、みちびき公式ページ内に公開されています。

QZS L6 Toolコード見直しの過程で、CLAS仕様書でわからなかったことを思い出しました。

  1. フレーム先頭の見つけ方。CLASが伝送されるL6D信号では、1秒間単位で1データパート（DP, 2000ビット）が伝送されます。受信機は、DPを5パート集めてサブフレーム（SF）を、さらに、SFを6サブフレーム集めて1周期をなすフレームを得て、その内容を解釈します。各DPのヘッダには、SF先頭を示すサブフレームインジケータ（1ビット長）があり、SF先頭は識別できます。しかし、フレーム先頭を識別する公式手段がないようです。QZS L6 Toolでは、SF先頭にマスク情報（サブタイプ1）があれば、そのSFをフレーム先頭とみなしています。マスク情報がないと、情報を復号できないからです。これが正しければ、仕様書4.1.2.2.1節（24ページ）にそのように書かれていると安心できます。
  2. 利用不可（not available）の扱い方。例えば、サブタイプ2のGNSS Hourly Epochの値が3600以上のときは利用不可を表すとされていますが、そのときも、続くSSR Update Intervalなどを読むことは自然だと思います。しかし、例えば、サブタイプ8のSTEC Polynomial CoefficientC00に利用不可を表すデータが格納されていたとき、続くC01やC10は無意味です。このときC01やC10にも利用不可を表すデータが格納されるのでしょうか。または、このときはC01やC10は無意味なので伝送しないのでしたら、現在のQZS L6 Toolではそれ以降のSFデータが読めなくなります。データ利用不可により無意味になる続きデータには、常に利用不可データが格納されることがNote部分に書かれていると嬉しいです。
  3. CLASでは、compact network IDとgrid IDから、仕様書4.1.4節（71ページ）にある、あらかじめハードコードされた参照座標値を求めます。はじめに、受信機は、この参照座標値群を用いて、観測したい場所の大まかな座標に最も近いnetwork IDとgrid IDを見つけ出すことになります。このとき、観測点座標値と参照座標値との距離のしきい値はどのように決定するのでしょうか。海上などでCLASを利用するとき、遠い参照点を利用しないようにするために、このしきい値の参考値が仕様書4.1.4節にあると嬉しいです。
  4. サブタイプ12のSTECデータの読み込み方法。STEC Correction Availability（STEC Correction Typeのことではありません）として2ビットが割り当てられています（仕様書4.1.2.2.13節、表4.1.2-35、52ページ）。しかし、仕様書56ページの(2)には、その1ビット目の動作しか記述されていません。この2ビット目は、reservedで、常にゼロがセットされるのでしたら、56ページの(2)にそのように書かれていると安心できます。
  5. STEC Quality IndicatorからSTEC Qualityへの変換方法。サブタイプ8でのこの説明として仕様書4.1.2.29節（43ページ）には「The definition is the same as that for SSR URA whereas the dimension is TECU instead of m」と書かれているので、URAの6ビット値から距離精度を求める式がそのまま利用できると解釈できます。ところが、5.4.2節の表5.4.2-1（81ページ）に書かれているURAの値と、5.4.3節の表5.4.3-1（84ページ）に書かれているTECUの値とが異なるのです。5.4.2節には値だけでなく計算式も書かれていますが、5.4.3節には計算式がありません。これらのURAを、ある値で割った値がSTECになる訳でもないようです。仕様書5.4.3節にSTEC Quality IndicatorからSTEC Qualityを求める計算式が記述されると嬉しいです。

CLASLIBのソースコードを読めば、これらの疑問点は解決できるかもしれません。CLASLIBは大きなソフトウェアでして、実は、私はそれを途中で読むのを諦めてしました。ソースコードが公開されていることはとてもありがたいことなので、少しづつ、読んでゆきたいと思っています。

CLASを利用すれば、その素晴らしさが実感できます。これが日本国内で無償で利用できることは、とてもありがたいことで、感謝しています。

まとめ

機械可読かつ人間可読を目指し、QZS L6 Tool出力を表形式にしました。また、その利用方法をまとめてみました。また、忘れないうちに、CLAS仕様書上の私の疑問点をまとめました。gitやgnuplotはとても便利です。

測位衛星は私たちの身近にあります。高精度測位補強サービスのCLAS、MADOCA-PPP、SLAS、Galileo HAS、BeiDou PPP-B2b、SBAS、は全部、素晴らしいです。その高度が技術が、誰でも無償で事前連絡なしに利用できることは、ありがたいことだと思っています。

関連記事