前回、「アルマナックデータによるGPS衛星の位置計算其の3」で3次元空間上でのアルマナックによるGPS衛星の位置(X,Y,Z)の求め方を掲載しました。今回は公約通り、これを経度緯度で表現する方法を一緒にまなんでいきましょう。これがわかれば、地球上のどのあたりの上空にGPS衛星があるかわかります。

尚、すでに、測地系の雑学をお話しした通り、いろいろなやり方がありますが、ここでは、一般的なWGS84に則した方法をとりたいと思います。（単に、小職がこの方法しか知らないというのも理由の一つですが・・・)

とはいいながら、軌道計算もさることながら、経度、緯度に変換する方法はさらにその上をいって難しいです。よってこれもまともに解くとわけが分からなくなります。（少なくとも小職はそうです)

というわけで、ここでは、WGS84測地系の定義を以下のように簡素化します。

○座標系

   ECEF系(ここは特に変更なしです。)

○準拠楕円体

楕円でなく真円とする。

○ジオイド面

高さは今回は意識しません。経度と緯度が求まればいいとします。

こんな感じです。実質、経度緯度を求めるという意味では、今回は楕円を真円と見なしているだけですが、実際は大問題です。しかしながら、だいたいこのへんと見るという前提で許してください。また、次ステップへの橋渡しと考えてください。

では、次から本格的に説明していきたいと思います。

[0回]

前回、アルマナックデータから3次元空間におけるGPS衛星の位置を求められましたので、そろそろ経度と緯度に関して触れておこうと思います。あえて、ふれてなかったところもありますが、経度と緯度はおくが深いです。

さて、次回あたりから、GPS衛星の位置を経度、緯度で表す方法をトライしてみようと思ってますが、その前に測地系について雑学を学んでおきましょう。

測地系という言葉は聞いたことはあるかとは思いますが、そもそも測地系って何？、正直、軌道計算を行う前までは、結構あいまいでした。定義自体は、

「地球上の位置を経緯度及び標高を用いる座標によって表すとき、前提とする条件」となっております。

つまり、測地系とは、経度、緯度を表現するために基とする条件なのです。

この測地系の条件として、3要素が上げられ、座標系、準拠楕円体、ジオイド面があります。

座標系は

中心を地球の重心を原点としたり、観測点を原点とするなど、座標の原点や、x、y、z軸の方向の定義に関わる部分です。

準拠楕円体は

地球は楕円ですが、その楕円の定義に関する部分です。

ジオイド面は

標高の基準の定義に関する部分です。

という感じですが、そもそも測地系って、なんで一つでないんだということになるのですが、今のように、技術が進歩してないときに、地図等をつくるために、いろいろな国が三角測量やら天体観測やらで、独自の測地系を作り上げたことが、理由にあるようです。

しかしながら、GPS等の発達によって、得られた情報による、正確な測地系の世界測地系へと変わりつつあります。

日本も、2002年に世界測地系に変わっています(思いのほか結構最近ですね)。とはいいつつも、例えば米国の世界測地系と日本の世界測地系とは違いがあります。

例えば、準拠楕円体は日本ではGRS80楕円体で、米国ではWGS84の楕円体だったりします。
米国の方だけ、「の」をいれましたが、このWGS84の中に、座標系やら準拠楕円体の定義が含まれてます。

しかしながら、ここで、一度つまづきました、なんで、日本はWGS84を採用しなかったのか、事実上、日本と米国の世界測地系にさほどの差異はないとのことですが、気になります。

ここら辺の回答というわけではないのですが、国土地理院のホームページにQ&Aがありそこを見るとなんとなくわかりました。

日本は世界測地系としてITRF系と呼ばれるものを採用しているそうなのですが、このQ&Aに採用した理由が以下のようになってました。

(1)　日本を含む国際協力で構築されている。我が国は、VLBI、GPS等の国際共同観測を通して、ITRF 系の構築に貢献しています。
(2)　高精度である。大量の宇宙測地観測結果から計算されるため高精度なものとなっています。
(3)　公開性が高い。一国により、又は軍事目的に構築されたものではないので、オープンなものとなっています。

なるほど、そういうことだったんですね。
政治的なところもあるようですね。

かなり大雑把でしが、ここらへんのネタはまたお話したいと思います。

[0回]

いよいよ大ずめです。前回「アルマナックデータによるGPS衛星の位置計算其の2」で離心近点角Eをニュートン法により算出しました。今回は、クライマックスです。

さてさて、ここまでくれば、後は気が楽です。

まず、前に「航法メッセージによる GPS衛星の位置の求め方 其の3」に記載した、平面軌道上のGPS衛星位置(x',y')を求める公式を利用して、2次元面での、GPS衛星位置を求めます。

これはもう簡単ですよね。

x' = a*cos(E) - a*e
y' = a*sqrt(1-e²)*sine(E)

これで算出できます。ここで、お分かりだとは思いますが念のため、aは長半径、eは離心率、Eは離心近点角になります。
あとは、アルマナックデータおよび、前回求めた離心近点角Eを当てはめて解くだけです。解くと以下のようになります。

x' = -11120449.437278792
y' = -24439569.638433773


さてさて、最後に、上記2次元座標を3次元座標に変換してやりますがこれも、「航法メッセージによる GPS衛星の位置の求め方 其の3」で記述した、回転行列を利用します。

まずは、行列を解きます。すると、

X     =  x' * cos(Ωk) * cos(ω) - x' * sin(Ωk) * sin(ω) * cos(i)
           - y' * cos(Ωk) * sin(ω) - y' * sin(Ωk) * cos(ω) * cos(i)

Y     =  x' * sin(Ωk) * cos(ω) + x' * cos(Ωk) * sin(ω) * cos(i)
      　　- y' * sin(Ωk) * sin(ω) + y' * cos(Ωk) * cos(ω) * cos(i)

Z     =  x' * sin(Ωk) * sin(i)　+ y' * cos(Ωk) * sin(i)


となります。少々見にくいですが、こんな感じになります。ここでΩkは前回求めた、現在時刻による昇交点経度、i、ωはどちらもアルマナックデータより取得できそれぞれ、軌道傾斜角、近地点引数になります。考え方のおさらいをしますと。

「3次元の軌道面は、2次元の軌道面をz 軸回りにω 回転させ，続いてx 軸まわりにi，z 軸回りにΩ 回転させたものである」

でしたよね。

さて、実際にときますと、以下のようになります。

X = -26312458.28576335
Y = -3442369.7618006254
Z = 4094091.9850355694

どうでしたか？とうとう3次元空間における、GPS衛星の位置情報が求まりました。これは紛れも無い解です。しかしながら、GPS受信機という意味では不十分です。といいますのは、航法メッセージの中にはアルマナックデータだけでなく、エフェメリスデータというものが入っており、こちらのほうが正確な情報を持っているのですが、残念ながら、このやり方では、エフェメリスデータに含まれる補正データがうまく扱えません。よって、新たな角度、真近点角を利用した方法に置き換える必要があります。え！じゃー最初からそれをおしえろよ！と言われるかもしれませんが、概念的には今までの内容を理解していれば、問題ありませんし、入門としては、今回の方法から入ったほうが分かりやすいと思います。

機会があれば、エフェメリスデータによるGPS衛星の3次元空間上での位置の求め方を公開したいと思います。

[0回]