2016年6月14日火曜日

GPS 経路データから求めた誤差について

六甲全山縦走路については,以前書いた六甲全山縦走路の距離についてや,六甲全山縦走路のデータや,後日書いた地理院地図で六甲全山縦走路の登山道が修正されていたも参考にして欲しい。

 これまで,GPS で取得した経路データを Google Map や OpenLayers を用いて地図上にプロットしたり,経路上で撮影したデジカメのデータを地図上に配置したりしてきた。 その際,GPS で取得した経路データは Douglas-Peucker のアルゴリズムで間引いている。そのデータの間引きについては以前このブログに記載した (GPSのログを間引く(折れ線を間引く))。 また,間引いた後の GPS 経路データから,ヒュベニの公式を用いて距離を求めたり,地理院地図のデータから各点の高度情報を取得してルートの高度情報を得たりしている。 また,移動距離と高度情報から,斜面に沿った沿面距離を求めたりもしている。

 この2年ほど,六甲全山縦走路を歩いている。 毎年,春に1回,秋に数回歩き,これまで個人的に8回完走している。 そこで気になっているのが六甲全山縦走路の距離の問題である。 六甲全山縦走路の距離についてもこのブログに記載したことがあるが,神戸市のサイトには 56 km とあるが, GPS 経路データから求めると,約 45 km にしかならない。 沿面距離にしてもたかだか 1 km 程度増えるか増えないか,である。 さすがに 45 km が 56 km にはならないので,56 km の根拠が不思議でならない。 以前,神戸新聞でも話題になったようであるが,六甲全山縦走路を歩く人たちは実質 45 km と思っているみたいである。

 それはさておき,ふと,GPS 経路データの誤差ってどれ位あるか,と思った。 そこで,ちょっとだけ考えてみた,というのが今回の話し。

 GPS 経路データは,生のデータは以下のような NMEA 形式のテキストファイルという形で得られる。 (これで1個の点を表している)
$GPGGA,054552.000,3404.6509,N,13611.5193,E,1,3,25.13,113.7,M,36.3,M,,*68
$GPGSA,A,2,13,04,17,,,,,,,,,,99.99,25.13,99.99*36
$GPGSV,2,1,07,04,62,340,34,17,62,162,24,13,48,100,33,20,08,053,*7F
$GPGSV,2,2,07,02,,,33,10,,,35,23,,,42*7C
$GPRMC,054552.000,A,3404.6509,N,13611.5193,E,12.17,357.20,261111,,,A*5C

 GPS だけでなく,ロシアのGLONASS の情報を受信できる GNSS ロガーのデータは以下のようになっている。
(これも1点のデータ。GPGSV は5回に1回の記録みたい)
$GPGGA,203323.000,3438.2725,N,13506.0287,E,1,11,0.98,19.1,M,33.5,M,,*50
$GNGSA,A,3,27,04,26,09,16,23,08,,,,,,1.75,0.98,1.45*1C
$GNGSA,A,3,75,85,73,74,,,,,,,,,1.75,0.98,1.45*16
$GPGSV,4,1,13,27,72,357,48,08,56,263,42,16,53,062,45,26,32,090,44*71
$GPGSV,4,2,13,07,24,317,,04,23,188,34,09,18,272,18,23,16,233,28*7A
$GPGSV,4,3,13,22,04,109,25,31,03,149,39,193,,,37,47,,,*4C
$GPGSV,4,4,13,01,,,30*79
$GLGSV,2,1,06,74,75,097,42,85,64,239,43,75,42,333,39,73,25,137,34*65
$GLGSV,2,2,06,86,12,221,28,83,06,034,*6F
$GNRMC,203323.000,A,3438.2725,N,13506.0287,E,0.18,69.09,231015,,,A*4A
$GPVTG,69.09,T,,M,0.18,N,0.34,K,A*05

 どちらにしても,GPGGA と GPRMC (GNRMC) という行の中に,緯度情報は「3438.2725,N」,経度情報は「13506.0287,E」という形で書かれている。 ここで,緯度情報の「3438.2725,N」だが,これは「北緯 34 度 38.2725 分」を表している。 これから緯度方向の距離の最小の刻み幅を見てみよう。

 0.0001 分が最小刻みとなっているが,角度で言うと 360 度で地球一周の約4万kmを表している。 すると,0.0001分 は,
0.0001/60/360 * 40000*1000 = 4.6e-9 * 4e7 = 0.185 [m]
に相当することになる。 つまり,我々が容易に手に入れることのできる GPS 経路データは 18.5 cm より細かくは刻んでくれない, 言い換えると約 20 cm の誤差を含むことになる。 (より精度の高い機械もあるようだが,それはきっともっと高価なはず)
 経度方向は,どの緯度にいるかによるが,一般には緯度方向よりも誤差は短い。 (cos(緯度)という因子をかけないといけない) しかし,日本ではせいぜい半分に減る程度なので,誤差を考える時にはあまり変わらない,という認識でいいと思う。

 実際に六甲全山縦走した際のデータを見てみると,六甲全山縦走路約 45 km に対して,GPS 経路点は約 44000 点である。 何度か記録させているが,毎回ほぼ同じ数となっている(それなりにバラつきはあるが…)。 私は1秒に1点記録させているので,45 km = 45000 m から,1秒に約 1 m 進んでいることを示している。 1秒に 1 m ということは,時速約 3.6 km という事を表している。 まぁ,そんなもんでしょ。

 このデータを間引くと,だいたい 2200 点前後になる。 つまり,隣り合う点間の距離は平均で約 20 m となる。 各点の位置精度が 20 cm の誤差を含むとなると,距離を求めると 1 % の誤差を含むことになる。 ここから言える事は,
六甲全山縦走路,約 45 km に対して,計測した距離データは最大で 400 m ~ 500 m の誤差を含んでいる
という事になる。 逆に言うと,誤差は 1 km は越えない,と考えられる。 この 400 ~ 500 m というのは最大であり,各点の誤差がプラスマイナスで打ち消し合うと,合計としてはもう少し誤差は減る可能性がある。 いずれにせよ,45 km という距離は誤差を考慮しても 56 km とは明らかに違う,ということを意味している。

 それでは,私が用いている距離計算はどの程度正しいだろうか?
比較のために,地理院地図上のツールと,国土地理院のサイトにある距離計算のページでの結果と比較してみよう。 念のため,南北方向と東西方向で比べてみた。 まず,比較した計算方法だが,以下の4つで比較している。
(a) 地理院地図の「機能」-「ツール」-「計測」
(b) 国土地理院測地部の下記サイト(距離と方位角の計算)
   http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2stf.html
   パラメーターは GRS80 を使用。
(c) 国土地理院測地部の下記サイト(距離と方向角の計算・ST計算)(100 km を越えると計算してくれない時があった)
   http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/xy2stf.html
   パラメーター:世界測地系,平面直角座標系=9系
(d) マイルーチン:
   http://hippo.matsup.mydns.jp/OLayerMap/ol3map_add_modify_line.html
 使い方はそれぞれトライして(調べてみて)欲しい。ちなみに,一番下のマイルーチンは,私が作ったもので,クリックで最初の点を指定し,次に終点をダブルクリックで指定する。 それから「Display Line Data」というボタンを押すと,下の方に2点間の距離を表示してくれる。
 上記サイトの距離の計算方法だが,(a) の地理院地図では「leaflet」という地図表示用の JavaScript の距離計算ルーチンを用いているらしい。 (b), (c) は各サイトに計算式が示されているので見てほしい。 (d) は私が作ったサイトであり,距離計算はヒュベニの公式を用いている。 (結局は国土地理院測地部から得た情報みたいやけど…)

 これらで2点間の距離を比較してみた。

(I) 銚子市役所前と日立駅前の距離
  まずは,南北方向の距離で調べたいので,銚子市役所前と日立駅前で比べてみた。 地図でそれぞれを拡大表示し,適当なポイントを決めて距離を求めている。 今回は銚子市役所前の道路の折れ曲がりの角と,日立駅前ロータリーにある噴水(と思われる)の中央の距離にしてみた。 なぜこの2点間かと言うと,おおまかに南北であり,距離が 100 km 未満の点だから,である。 別に他の2点間でもよかったが,最初に銚子市役所から,と考えたので,この2点間にしてみた。 100 km 未満,というのは (c) のルーチンが長い距離だと計算してくれなかった,からである。
 以下がその結果である。
 (a) 地理院地図:96.649 km (マイ/地理院地図 = 0.99689598): -0.31 %
 (b) 国土地理院測地部:96.3487 km (マイ/地理院 = 1.00000311,δ=3.1e-6)
 (c) 国土地理院測地部ST:96.3493 km (マイ/地理院ST = 0.999999969,δ=-0.31e-7)
 (d) マイルーチン:96.349 km
 結果を見ると,いずれもほぼ一致している。特に (b), (c), (d) の間の一致がとてもよい。 それに比べると leaflet のルーチンを用いている (a) だけ多少異なってる。 それでも 0.3 % 程度の違いなので,GPS 経路データを用いる上では問題ないと思われる。

(II) 銚子市役所前と東京駅前の距離
  次に同じ銚子市役所前から東京駅前の距離を求めてみた。これはほぼ東西の方向の距離の比較になっている。
 (a) 地理院地図:96.150 km (マイ/地理院地図 = 1.0011128): 0.11 %
 (b) 国土地理院測地部:96.2571 km (マイ/地理院 = 0.999998961,δ=-1.0e-6)
 (c) 国土地理院測地部ST:96.2574 km (マイ/地理院ST = 0.99999584,δ=-4.1e-6)
 (d) マイルーチン:96.257 km
 これもよい一致をみている。特に下の3つは excellent といいたくなる一致である。 一方 (a) はやはり多少ずれているが,それでも 0.1 % 程度なので,良い一致と言える。

 では,もっと距離が遠くなるとどうだろうか?
実は,当初は銚子と稚内,銚子と対馬,で計算させたのだが,手法 (c) は,距離が長いと計算してくれなかったので,上記の2パターンにしてみた。 折角なので,銚子と稚内,銚子と対馬,の結果も載せておこう。

(III) 銚子市役所と稚内港の突堤の距離
  以下は銚子市役所前の道路の角と,稚内港の開運二丁目の突堤の右上角との距離である。 (I),(II) よりも遠方で,約 1000 km 離れた2点間である。
 (a) 地理院地図:1080.120 km (マイ/地理院地図 = 0.9975725) :- 0.25 %
 (b) 国土地理院測地部:1077.4854 km (マイ/地理院 = 1.0000117):0.001%
 (c) 国土地理院測地部ST:距離過大で不可
 (d) マイルーチン:1077.498 km
これより,マイルーチン(ヒュベニの公式)と国土地理院のサイトの差はほぼ無視できる。 一方,leaflet のルーチンを用いた地理院地図では,差が 0.25 % となっている。 これは 100 km 程度の場合と大差ない結果となっている。

(IV) 銚子市役所と対馬の北の三ッ島灯台との距離
  こちらは銚子市役所前と対馬の少し北にある三ッ島にある灯台の間の距離である。 こちらも約 1000 km 離れているが,東西方向となっている。
  地理院地図:1040.523 km (マイ/地理院地図 = 1.001639) :0.16 %
  国土地理院測地部:1041.6321 km (マイ/地理院 = 1.000572):0.06 %
  国土地理院測地部ST:距離過大で不可
  マイルーチン:1042.228 km
結果は,国土地理院測地部との差は増えているが,それでも地理院地図との差 0.1 % よりは小さい。 しかし,距離が長くなると,どうもヒュベニの公式だとズレが大きくなるみたいな気がするなぁ。

 上記の4パターンの比較から,私が距離計算に用いているヒュベニの公式は,国土地理院測地部の複雑な計算方法とまずまずの一致を示していることがわかった。 但し,距離が 100 km を越える時には多少誤差が増えるみたいだった。 一方,同じ国土地理院が提供している地理院地図のツールでは,0.1~0.3 % 程度の誤差があることがわかった。 leaflet のルーチンはその程度の誤差が含まれる,ということみたい。

 上記のことから,私が GPS 経路データから求めた六甲全山縦走路の距離は,500 m 以内の誤差で約 45 km となる。 念の為に GPS ロガーのデータではなく,手で地図上に点をプロットして距離を求めてもほぼ 45 km だった。 やはり,神戸市はなぜ公式サイトに 56 km と書いているのか,その根拠を示して欲しい,と思う,今日この頃である。

0 件のコメント: