[Peer Reviewed]
Proceedings of the 5th Conference of the Society for the Advancement of Positioning Technologies (SAPT), Osaka, Japan, Sep. 1, 2023

GNSS 解析ソフトウェアから求めた対流圏遅延量推定の比較
Comparison and evaluation of tropospheric delay estimation
obtained from GNSS analysis software
中川 豊 † ，東野 武史 † , 岡田 実 † ,† 奈良先端科学技術大学院大学 先端科学技術研究科
Yutaka Nakagawa† , Takeshi Higashino† and Minoru Okada†† Nara Institute of Science and Technology
Graduate School of Advanced Science and Technology
Abstract The propagation delay caused by radio waves from GNSS satellites passing throughthe troposphere is treated as a tropospheric delay and an error factor during positioning. On
the other hand, if the coordinates of observation points are accurately measured, the amount of
delay due to the troposphere can be accurately estimated, and it is applied to weather prediction
as an indicator of water vapor in the atmosphere. Tropospheric delay is also expected from theperspective of disaster prevention in response to heavy rain disasters, which have become more
frequent in recent years.GNSS analysis software requires more accurate satellite positioning, butin this paper we focused on tropospheric delay and veriﬁed the performance of tropospheric delay
obtained from multiple academic GNSS analysis software.
キーワード: GNSS, 対流圏遅延, RTKLIB, G-RitZ

1

はじめに上空の水蒸気の急増は大雨の発生前の予兆のサインと見なすことができ, 水蒸気の正確な把握により線状降水帯のような大雨の予測精度の向上が期待されている. 局地的な気象予測には空間分解能の向上が前提であり,衛星測位の誤差要因の１つである対流圏遅延の高精度かつ迅速な推定が重要であると考えられる. 対流圏遅延量は大気の水蒸気の状態を表現しており, 天気の変化に対応している. しかしなから現時点では大気中の水蒸気を精度よく, リアルタイムかつ広域 (国土地理院の GEONET を利用) に把握できる方法は GNSS 測位情報から推測される対流圏遅延量以外にはない. 対流圏遅延量は GNSS では直接観測できないため精度のよいモデル化, 推定が必要である. さらに対流圏遅延量の推定精度の向上は測位精度の向上につながる. 本稿では GNSS測位解析ソフトウェアの違いは, 対流圏遅延の推定に対してどの程度精度に差があるのかに注目し解析, 検証を行った.

2

各解析ソフトウェアの測位手法
本稿執筆時 (2023 年 10 月 13 日現在) において, 学術調査研究用として世界的な定評と実績のある主な測位解析ソフトウェアは Bernese Version5.2( University of Bern)[1], GAMIT/GLOBK Release10.4(Massachusetts Institute of Technology) [2], GIPSY-OASIS II(Propulsion Laboratory)[3] および RTKLIB
v2.4.3 (Tokyo University of Marine Science and Technology) [4] などがある. 測位手法は, スタティック測位／キネマティック測位, ネットワーク測位／短基線／ PPP（Precise Point Positioning）に分類され, 主な解析ソフトウェアの測位手法の対応を示したのが表 1 である.

第5回測位技術振興会研究発表講演会
講演論文集（2023年9月1日，大阪）

-21-

表1

各解析ソフトウェアの測位手法の対応

GAMIT

GIPSY-GLOBK-OASYS〇〇〇〇〇〇〇〇〇RTKLIB

G-RitZ〇〇※〇〇〇〇※〇〇〇〇〇※〇〇〇〇〇※〇

Bernese
Networked
static positioning
Networked
kinematic positioning
Single- baseline static
relative positioning
Single-baseline kinematic
relative positioning
PPP
Static positioning
PPP
Kinematic positioning※ Possible (neary) Real-time processing3

大気遅延誤差とその推定モデル

3.1

大気遅延の原理

大気屈折による大気遅延は２つの物理的効果に起因する. １つは, マイクロ波が誘電媒質である大気中を通過するために真空中よりも減速され, 見かけ上伝搬経路が伸びる効果である. ２つめは, マイクロ波の伝搬経路が曲率 (raybending) を持つため, 直線より経路が実際に長くなる効果である. これらの効果による遅延 ∆L は次式で表される [5].

∫[n (s) − 1] ds + [S − G]

∆L =(1)

L

ここで,G は GNSS 衛星から地上のアンテナまでの直線距離,n(s) は大気屈折によって曲率が生じた伝搬経路

L 上の点 s における屈折率である (図 1). このとき,S は L 上の微少部分 ds を積分して得られる. 右辺第１項が減速の効果であり, 実際の電波は L 上を大気中の伝搬速度 catm で伝搬するが, 解析の上では真空中の速度 c で伝搬し, 見かけ上伝搬経路が伸びたことになる. そして, 残りの右辺第２項が伝搬経路の曲率の効果となり, この効果による遅延は天頂方向で 0 となる. 衛星 p の電波信号を受信機 u で観測して得られる距離測定値は次のように示される.

f12 p
δI (t) + δTup (f ) + epi,u (t)fi2 u(2)f12 ppδI (t) + δTup (t) + λi Ni,u+ εpi,u (t)fi2 u(3)ρpi,u (t) = rup (t, t − τup ) + c {δtu (t) − δtp (t)} +ϕpi,u (t) = rup (t, t − τup ) + c {δtu (t) − δtp (t)} −こ こ で 式 (2) は コ ー ド 擬 似 距 離 [m], 式 (3) は 搬 送 波 位 相 観 測 量 [m] で あ る. た だ し,i,fi , λi は そ れ ぞれ 搬 送 波 の 周 波 数 帯 を 表 す 番 号, 周 波 数 お よ び 波 長 で あ る.

GPS で は, i = 1, 2, 5 に 応 じ て, fi =1575.42, 1227.60, 1176.45 [MHz] が割り当てられている. このうち δTup が本稿で議論する, 衛星 p と受信機アンテナ u の間に存在する対流圏遅延量を表す. ここで rup は衛星 p と受信機 u の幾何学的距離であり,p − u 間の電波伝搬距離を τup , 光速を c, 受信機時計誤差を δtu , δtp は衛星時計誤差を, δIup は L1 帯のコード
ppp

擬似距離に対する電離圏遅延量を, Ni,u は整数値バイアスを, ei,u , εi,u は観測雑音を表す.

-22-図13.2

図2大気屈折の模式図

マッピング関数の模式図マッピング関数

正確な大気遅延を得るためには, 衛星と受信機のアンテナ間の伝搬経路上の気温・気圧・水蒸気分圧の値から大気屈折率を計算し, これを積分する必要がある. しかしながら大気遅延は時間・空間の双方の関数でありGNSS の解析では, このままで大気遅延の観測方程式を解くことは事実上不可能となる. そこである一定の時間内において大気が水平方向にはほぼ一様な構造を維持するとすれば, 遅延量は仰角のみに依存する関数となり, 未知数を大幅に減らすことができる. この仰角依存関数をマッピング関数 (mapping function)と呼び衛星と受信機のアンテナ間のある仰角 esr の全遅延量 Trs は下記のように表すことができる (図 2)

Trs = Mh (esr )ZHD + Mw (esr )ZWD

(4)

ここで, ZHD と ZWD はそれぞれ天頂方向の 静水圧遅延量 , 湿潤遅延量 であり, ZTD = ZHD + ZWDの関係にある. Mh と Mw がそれぞれのマッピング関数であり次式のように表される [6].

Mh =1sin (esr ) + tan(e0.00143s )+0.00445

(5)r

Mw =

1(6)

0.00035sin (esr ) + tan(es )+0.017
r[7] によると,[6] は対流圏遅延量の補正誤差が日本付近で大きいことや, マッピング関数の性能は湿潤大気の精度に依存するため，元となる数値気象データの水蒸気分布の正確な把握の必要性などが指摘されているが, 本研究では一般に広く知られている数値予報データから構築された簡便な経験的マッピング関数である [6] のGMF（Global Mapping Function）モデルのマッピング関数を採用した. なお [7] の指摘した内容については,今後再評価を検討している.

3.3

GNSS による対流圏遅延と可降水量の計算

GNSS から得られた ZTD から, 式 (7) によって地上気圧から求められる ZHD を差し引いた ZWD を式 (8)に当てはめると GNSS 由来の天頂可降水量 (以降 GNSS-pwv) が求められる.

ZHD =k 1 · Rd· P0 × 10−6g(P W V = ZW D ×Rv

105
k3k2 +Tm(7))(8)

ここで k1 , k2 , k2 はそれぞれ 77.604 [K/hPa], 64.79 [K/hPa], 3.776×105 [K2 /hPa], Rd は乾燥空気の気体定数, Rv は水蒸気の気体定数, g [m/s2 ] は重力加速度, P0 [hPa] は地上気圧である. 大気の気温は直接観測す-23-

る方法の他、地上気温 Ts [◦ C] を用いた式 (9) を利用して大気の加重平均気温 Tm として代用する.

Tm = 70.2 + 0.72(Ts + 273.15)(9)

なお Tm は [8] に示すおり経験式であり，ラジオゾンデの観測データから求めた水蒸気分圧で重み付けられた気柱の加重平均された気温を示す.

4

測位解析ソフトウェアの違いによる対流圏遅延の推定の比較
表2

解析ソフトウェアのパラメータ設定

Program version

RTKPOST ver.2.4.3 b34
G-RitZ ver.0.1-bata

Obs data

GEONET YONAGO A station

Nav data

GEONET YONAGO A station

Ephemeris data

IGS sp3

Clock data

gbm clk

Ion data

not useObs start

2017/11/1 00:00 UTCObs end

2017/11/30 23:59 UTCFrequency

L1+L2

SatellitesGPS

Est method

Kalman ﬁlter

Pos mode

PPP StaticSolution

ForwardElev mask

15.0 deg

Weighting of obs dataoﬀ

by Elev maskDynamicsoﬀ

Tide corron

Ion. corr

Ion free lc

Trop. corr

Saastamoinen model

Trop. opt

Estimate ZTD

Phase Center Variation

GSI PCVDiﬀ. Code Biason(Receiving)

対流圏遅延による降雨予測を行う場合, 対流圏遅延推定の精度が重要と考えられるが, 測位解析ソフトウェアが及ぼす対流圏遅延推定値の違いや誤差に関する評価は不十分である。そこで本稿では, 国内を中心に広く利用されているオープンソースの高精度測位用プログラムである RTKLIB [4] と, G-RitZ [9][10] が導出する対流圏遅延を比較する. 表 2 は RTKLIB の後処理ライブラリである RTKPOST と G-RitZ の対流圏遅延推定の精度を比較するために設定した主なパラメータを示す. 図 3 および図 4 は, それぞれ表 2 に基づき 2017年 11 月 2 日の 1 日間における,RTKLIB と G-RitZ の ZTD,ZHD および ZWD の時系列である. 一方, 気象庁は上空の気圧, 気温, 湿度, 風向, 風速等の気象要素を観測するラジオゾンデという気象観測測器による観測されたデータがある [11]. 得られた観測値を基に地上から上空の水蒸気の積分として可降水量 (sonde-pwv) を求めることができる. 直接水蒸気を観測できる水蒸気ライダーは, 高精度であるが, 現状まだ試験的運用の段階-24-である. 当日の気象状況として米子に近い松江の高層気象観測データと対流圏遅延推定の関係を対比してみた.
図 3 と図 4 を比較してみると,ZHD は両者ともにほとんど変化が見られない. ZHD は式 (2) のように地上の気圧値が与えられれば，高い精度で計算することができ, ZHD は ZTD の中で大きな割合を占めるが，その変動量は小さく, かつゆるやかである特性を示している．一方で ZWD は数 cm〜数十 cm 程度であるが，場所や季節，天候によって大きく変化することが特徴的である. 図 5 は 2017 年 11 月 1 日から同月 3 日にかけての松江の上空の高層気象観測データを断面図で表したものである (気象庁提供). 横軸は時間軸で 00/02 は同月

2 日 00 時,12/02 は同日 2 日 12 時 (UTC) である. 縦軸は上空の高度を気圧を用いて表現したもので (例えば850hPa は上空約 1500m 付近), 色の分布は QV という水蒸気の混合比 (湿潤空気に含まれる水蒸気質量と湿潤空気質量の比) を表す. 図 5 の赤色の折れ線グラフは TQV であり, sonde-pwv を示す.TQV は地表から大2気上端までの単位断面積の気柱に含まれる水蒸気の総質量 [kg/m ] で, PWV[mm] と同義である. 00/02 から12/02 にかけて 830hPa から地上付近にかけて水蒸気の割合が高い部分が表現されている. 図 5 と図 3,4 を見比べてみると,G-RitZ の方が ZWD と水蒸気の変化の対応がよく表現できている.

図 3 RTKLIB

図 4 G-RitZ

図 5 sonde20171102-25-5

可降水量との比較
2017 年 11 月の 1 ヶ月間におけるラジオゾンデの観測によって求められた可降水量 sonde-pwv と GNSS 由来の可降水量 (GNSS-pwv) を比較した (図 6,7).
散布図 6,7 のそれぞれ横軸は sonde-pwv, 縦軸は GNSS-pwv である. 図中の直線は線形回帰直線である. ラジオゾンデは 2 回/日,00,12(UTC) に観測が行われるが, 観測時刻の 30 分前にラジオゾンデが放球されることから GNSS-pwv も 00,12(UTC) の前 30 分間の平均値を使用した. 理想的 (決定係数 r 2 = 1 の場合) には散布図としては, y = x の線上付近にプロットが集中するはずであるが, 両者ともにバラツキが大きく, 特にRTKLIB は sonde-pwv と GNSS-pwv の回帰直線の傾きがほとんど 0 という結果となった.

図 6 RTKLIB-pwv

図 7 G-RitZ-pwv

6

まとめと今後の課題
受信機や解析ソフトウェアによる測位精度の違いについて論じた研究は数多くみられるが, 本稿のように対流圏遅延および可降水量を比較, 検証した例はあまりない. 本稿では 2 つの測位解析ソフトウェアを用いて対流圏遅延および可降水量の比較を行った. その結果, ソフトウェアにより対流圏遅延および可降水量の精度に差が見られた. これはカルマンフィルタの適用方法の違いによるものと推測される. RTKLIB 並びに G-RitZ-26-

のカルマンフィルタに設定したシステムノイズの標準偏差は, ともに 10−4 とし, 両者に差を設けなかったが,事前誤差共分散行列の省略時標準の初期値は, 前者は (0.3, 0.3, 0.3), 後者は (0.3, 0.03, 0.03) であり予測ステップにおける違いがあるため, 回帰精度の差の一因となったと考えられる.
なお本稿では対象とした観測点は 1 か所であり, 観測期間も 1 ヶ月間であった. 今後の課題として季節的, 空間的に更なるデータを蓄積し, 検証の精度を高める必要があると考えている.

7

謝辞
本研究において立命館大学情報通信システム研究室から G-RitZ の提供をいただきここに感謝の意を表します.

参考文献[1] Rolf Dach and Peter Walser. Bernese gnss software version 5.2. 2015.
[2] TA Herring, RW King, SC McClusky, et al. Introduction to gamit/globk. Massachusetts Instituteof Technology, Cambridge, Massachusetts, 2010.
[3] Adem G Hayal and D Ugur Sanli. Revisiting the role of observation session duration on precisepoint positioning accuracy using gipsy/oasis ii software. Boletim de Ciências Geodésicas, Vol. 22,pp. 405–419, 2016.
[4] Rtklib: An open source program package for gnss positioning. https://rtklib.com/. (Accessedon 01/10/2024).
[5] 市川隆一, 関戸衛, 小山泰弘. 1d15 数値天気予報データを用いた大気伝搬遅延量推定ツールの開発. 2006.
[6] Arthur E Niell. Global mapping functions for the atmosphere delay at radio wavelengths. Journalof geophysical research: solid earth, Vol. 101, No. B2, pp. 3227–3246, 1996.
[7] H Munekane, Y Kuroishi, Y Hatanaka, and H Yarai. Spurious annual vertical deformations overjapan due to mismodelling of tropospheric delays. Geophysical Journal International, Vol. 175,
No. 3, pp. 831–836, 2008.
[8] Michael Bevis, Steven Businger, Steven Chiswell, Thomas A Herring, Richard A Anthes, ChristianRocken, and Randolph H Ware. Gps meteorology: Mapping zenith wet delays onto precipitablewater. Journal of Applied Meteorology (1988-2005), pp. 379–386, 1994.
[9] Atsushi Mouri, Yoshifumi Karatsu, Goshi Okuda, Sueo Sugimoto, Yukihiro Kubo, and Masaharu
Ohashi. New ppp/vppp algorithms by using multiple antennas. Transactions of the Institute ofSystems, Control and Information Engineers, Vol. 29, No. 12, pp. 525–534, 2016.
[10] Gnss 測位ソフトウェア g-ritz – 測位技術振興会. https://jsapt.net/ja/post/17070006. (Accessedon 01/10/2024).
[11] ラジオゾンデによる高層気象観測. https://www.jma.go.jp/jma/kishou/know/upper/kaisetsu.
html.

-27-