s-taka.org

高橋 賢

Lang:

準天頂衛星みちびきの測位ライブラリCLASLIBの利用（その2：ssr2osr）

category: gnss
tags: clas qzss

状態空間表現（SSR）

前回のCLASLIB構成の理解の続きです。今回は、CLASLIBのssr2osrを試します。

SSRとは、状態空間表現（state space representation）のことです。これは、衛星測位に関する状態を、例えば衛星内部時計誤差をより正確に表現することにより、測位精度を高める方法です。今回の精密単独測位（PPP: precise point positioning）はこのSSRに分類されます。

一方、OSRは、観測空間表現（observation space representation）のことであり、信頼できる別の観測所での観測データに対する補正値を利用者が得ることにより、利用者側での測位精度を高める方法です。かつてのDGPS（differential GPS）や、近年、期待されているRTK（realtime kinematic）は、このOSRに分類されます。

SSRを実現する具体的なプログラムとして、ssr2osr（SSR-to-OSR conversion）があります。CLASLIBでのssr2osrは、利用者側において利用するソフトウェアであり、CLAS信号のSSR補正により利用者の観測データをOSRに変換します。そのOSRを用いて、別のソフトウェアにて測位計算を行います。Geo++社Dr. Martin Schmitzのスライドの11ページにあるように、ssr2osrは、補正情報作成側で実行しても、利用者側で実施しても良いのです。最終的には、利用者側が、ssr2osrを実行することなく、このSSRを直接的に解釈することが望ましいとされています。CLASLIBでのrnx2rtkpは、このSSRの直接解釈を行うものだと思います。

自らの観測データを使ったCLASLIB ssr2osrの実行

ここでは、CLASLIBにて実現されるssr2osrを使います。

はじめに、CLASLIBをダウンロードして、そのutilsディレクトリ下のssr2osrディレクトリにてmakeを実行して、アプリケーションプログラムssr2osrを作成します。以下の作業は、このssr2osrディレクトリで行います。私は、LinuxのDebian 10（buster）ディストリビューションを利用しました。

次に、RINEX（ライネックス）ファイルである、衛星信号の観測擬似距離（OBS: observation）データファイルと、衛星軌道（NAV: navigation）データファイルを用意します。これらを国土地理院の電子基準点データ提供サービスから取得しても良いのですが、今回は私のRTK基準局データを利用しました。NTRIP Casterからインターネット経由でこのRTK基準局のRTCM3ファイルを取得し、そのRTCM3ファイルからRINEXファイルに変換します。このデータ取得と変換のためにRTKLIBを使います。

RTKLIBのソースコードをダウンロードして展開します。そして、そのappディレクトリでmakeを実行して、コマンドラインツールstr2strを作成します。このstr2strをこのssr2osrディレクトリにコピーした上で、

./str2str \
  -in ntrip://ntrip.phys.info.hiroshima-cu.ac.jp:80/F9P \
  -out file://%Y%m%d%H.rtcm3::S=1 -f 30

を実行すると、カレントディレクトリに年月日と時刻コードを表すファイルが作成されます。拡張子はrtcm3です。この時刻コードは、0時から23時までの時刻をaからxまでの英子文字で表したものです。str2strのオプションS=1 -f 30は、前後30秒間の余裕を持って1時間ごとに別のファイルを作成するものです。ここでは、RTK基準局で観測した2020年5月5日の3時から4時までのデータ20200505d.rtcm3（約3.8 MB）を用います。ここでの時刻表現は、UTC（世界標準時）でして、JST（日本標準時）よりも9時間だけ遅れています。したがって、このデータは、JSTでは12時から13時までのものです。

次に、このRTCM3ファイルからOBSファイルとNAVファイルを作成します。この変換にもRTKLIBにあるconvbinを利用します。

./convbin 20200505d.rtcm3 -v 3.02 -os -o 20200505d.obs -n 20200505d.nav

オプション-v 3.02はRINEXファイルバージョン3.02形式で出力することを（これを指定しないとバージョン2.11形式になります）、また-osは信号対雑音電力比（SNR: signal-to-noise power ratio）を結果に含めることを表します。CLASLIB利用時には、電波強度が一定以下の衛星信号を利用しないように設定するため、RINEXファイル作成時にこのオプション指定が必要です。オプション-oは、出力OBSファイル名指定のための、また、-nはNAVファイル名指定のためのものです。

なお、RTKLIBは、このような変換後のgzip圧縮ファイル（例えば、20200505d.obs.gzや20200505d.nav.gz）も、そのまま読み込めます（2021-03-19追記）。

次に、その日時に対応するCLASデータをダウンロードします。CLASアーカイブにて、UTCでの日時を選択して、ダウンロードします。CLASアーカイブでのファイル名は、年、1月1日からの日数、および時間コード（ただし、先の小文字を大文字で表す）からなります。この日付2020年5月5日は1月1日から数えて126日目なので、そのファイル名は2020126D.l6です。アーカイブページから日付を選択してファイル名をクリックしてそのファイルを取得する代わりに、cURLコマンド

curl https://sys.qzss.go.jp/archives/l6/2020/2020126D.l6 -o 2020126D.l6

にてそのファイルを取得しても良いでしょう。

そして、次のようなssr2obs設定ファイルosr.confを用意します。

元のサンプル設定ファイルからの変更点は、(1) 受信アンテナ形式指定をしないこと、(2) このRTK基準局の受信周波数帯に合わせて2周波数（L1帯とL2帯）にしたこと、(3) 出力ファイルをNMEA形式でなく緯度経度楕円体高（LLH: latitude longitude height）形式にして結果をわかりやすく表示したこと、の3点です。

そして、ssr2osrを実行します。

./ssr2osr -k osr.conf 20200505d.obs 20200505d.nav 2020126D.l6 -o 20200505d.llh

結果

作成された出力ファイル20200505d.llhの内容は次の通りです。このファイルには多くの情報が並ぶので、その最初の部分のみを抜粋します。

% program   : RTKLIB ver.2.4.2
% inp file  : 20200505d.obs
% inp file  : 20200505d.nav
% inp file  : 2020126D.l6
% obs start : 2020/05/05 02:59:12.0 UTC (week2104 183570.0s)
% obs end   : 2020/05/05 04:00:11.0 UTC (week2104 187229.0s)
% pos mode  : ssr2osr
% solution  : forward
% elev mask : 15.0 deg
% dynamics  : off
% tidecorr  : on
% tropo opt : off
% ephemeris : broadcast+ssr apc
% navi sys  : gps galileo qzss
% antenna1  :                       ( 0.0000  0.0000  0.0000)
%
% (lat/lon/height=WGS84/ellipsoidal,Q=1:fix,2:float,3:sbas,4:dgps,5:single,6:ppp
,ns=# of satellites)
%  UTC                   latitude(deg) longitude(deg)  height(m)   Q  ns   sdn(m
)   sde(m)   sdu(m)  sdne(m)  sdeu(m)  sdun(m) age(s)  ratio
2020/05/05 02:59:12.000   34.440129246  132.414755496   233.7646   5   9   4.776
4   3.5026   7.3992  -1.7312  -3.5717   3.4371 10000.00    0.0
2020/05/05 02:59:13.000   34.440130299  132.414757135   233.3983   5   9   4.775
8   3.5023   7.3982  -1.7303  -3.5711   3.4354 10000.00    0.0
... snip ...

これは、UTC時刻、緯度、経度、楕円体高、品質(5はPPP単独測位）をそれぞれ表します。品質が6（PPP）にならない理由はわかりません。(2022-04-01更新)

その第14列目のage(s)は、補強データが届いてからの経過秒数が記述されています。CLAS補強情報は5秒ごとに送信されますので、この結果の14列目の値が0, 1, 2, …5の範囲であるならば、CLASによる継続的なPPP測位が行われたことになります。この列の値が最大値の1000であれば、CLASによる測位の補強がなされなかったことになります。この結果から、時刻、緯度、経度、楕円体高、経過秒数のみを抽出するために、awkを用います。

awk '{print $2 "\t" $3 "\t" $4 "\t" $5 "\t" $14;}' 20200505d.llh

このCLAS情報ファイル2020126D.l6は、UTC時刻03:00:00から04:00:00までのものです。この時刻範囲で補強されれば、CLAS-PPP測位ができたといえます。

UTC     latitude(deg)   longitude(deg)  height(m)       age(s)
02:59:12.000    34.440129246    132.414755496   233.7646        10000.00
02:59:13.000    34.440130299    132.414757135   233.3983        10000.00
... snip ...
02:59:57.000    34.440126204    132.414755699   235.2483        10000.00
02:59:58.000    34.440127141    132.414755420   235.3133        1.00
02:59:59.000    34.440128971    132.414755250   235.2460        2.00
03:00:00.000    34.440131533    132.414754938   235.1163        3.00
03:00:01.000    34.440133495    132.414752663   235.4649        4.00
03:00:02.000    34.440134454    132.414753118   235.1127        5.00
03:00:03.000    34.440133393    132.414753927   234.9894        1.00
03:00:04.000    34.440134986    132.414753337   234.8666        2.00
03:00:05.000    34.440134555    132.414752899   234.7548        3.00
... snip ...
03:59:31.000    34.440113440    132.414782134   234.0389        4.00
03:59:32.000    34.440113328    132.414782015   234.2112        0.00
03:59:33.000    34.440113661    132.414782142   234.2457        1.00
03:59:34.000    34.440115799    132.414783428   234.1545        2.00
03:59:35.000    34.440116548    132.414784360   234.0657        3.00
03:59:36.000    34.440115588    132.414785932   234.2324        4.00
03:59:37.000    34.440114654    132.414787330   234.3477        0.00
03:59:38.000    34.440112902    132.414788481   234.3745        1.00
03:59:39.000    34.440112801    132.414785697   234.0207        2.00
03:59:40.000    34.440112910    132.414783535   233.6221        3.00
03:59:41.000    34.440113235    132.414783618   233.7450        4.00
03:59:42.000    34.440112895    132.414782548   234.0283        5.00
03:59:43.000    34.440112306    132.414782202   233.9166        6.00
03:59:44.000    34.440111882    132.414781813   233.8585        7.00
03:59:45.000    34.440110743    132.414781648   233.8334        8.00
03:59:46.000    34.440112951    132.414781696   233.9538        9.00
03:59:47.000    34.440111518    132.414779590   233.8481        10.00
03:59:48.000    34.440110554    132.414779131   233.9800        11.00
03:59:49.000    34.440110254    132.414777619   234.0704        12.00
03:59:50.000    34.440113086    132.414777736   233.5080        13.00
... snip ...

結果は、予想よりも2秒早い02:59:58から補強開始され、予想よりも24秒早い03:59:37の補強を最後に補強信号が利用されていないことを示しています。UTC時刻03:59:37以降は、その時刻での補強情報が利用され続けていることを表しています。

経度1度の差は約111キロメートルに相当します。静止しているRTK基準局のデータを処理した結果は、緯度の小数点以下5桁目まで固定されていますので、数十センチの精度と見積もることができます。この結果はあまり良いものではありませんが、このCLAS測位方法がPPP-RTKではなくPPPなので、そのようなものなのかもしれません。

これらの手順をシェルスクリプトにしました。下記をclas.shなどの名称でssr2osrディレクトリに保存し、実行属性をつけます（chmod 755 ./clas.sh）。また、このディレクトリにRTKLIBのconvbinを置きます。このclas.shを実行すると、上述の結果が得られるはずです。

#!/bin/bash
CLASF=2020126D
GNSSF=20200505d
# configuration file download
curl -s https://gist.githubusercontent.com/yoronneko/430ad5519e0d30a503eb643b349b8796/raw/4e5b7a0cbb1c703046e9ca97bd210c2bcb92b105/osr.conf -o osr.conf
# RTCM3 file download
curl -s https://s-taka.org/img/$GNSSF.rtcm3 -o $GNSSF.rtcm3
# CLAS augmentation file download
curl -s https://sys.qzss.go.jp/archives/l6/2020/$CLASF.l6 -o $CLASF.l6
# extraction of the observation and navigation files from the RTCM file
./convbin $GNSSF.rtcm3 -v 3.02 -os -o $GNSSF.obs -n $GNSSF.nav
# execution of CLAS PPP positioning
./ssr2osr -k osr.conf $GNSSF.obs $GNSSF.nav $CLASF.l6 -o $GNSSF.llh
awk '{print $2 "\t" $3 "\t" $4 "\t" $5 "\t" $14;}' $GNSSF.llh

次回は、CLASLIBを用いたVRS（仮想基準局: virtual reference station）によるRTK測位を報告します。

関連記事