先日、WGS84測地系についていろいろ掲載したが、たいてい、緯度経度について、勉強し始めると、日本で暮らしている人ならば日本測地系という言葉を耳にする。

測地系にはいろいろ、国によって独自の進化をとげたものがあるが、日本測地系もその一つです。

この測地系は衛星も存在しない、明治時代に形ができていたようで、準拠楕円体はベッセル楕円体というものを採用し、緯度経度は、天文観測で求めていたというある意味神秘的な話です。

ベッセル楕円体は赤道半径を6377397.155mとし、逆扁平率を299.1528128と定めているもので、ドイツのフリードリヒ・ヴィルヘルム・ベッセル(1784-1846)が計算して、その後、ドイツと日本などがこの基準を採用したとのこと。

前にも掲載しましが、日本は2002年までこの測地系でやってましたが、衛星や技術の進歩により、より正確な測地系に変わったわけです。

尚、世界測地系の経度緯度にあわせるならば、北西方向へ約450mずれるようです。

[0回]

前回の「WGS84測地系モドキによるGPS衛星の経度緯度の求め方其の2」の続きです。前回、座標系の変換を行いましたので、今回は、地球を楕円でなく球とみなして、いよいよ経度と緯度を求めてみます。

正直、球とみなせば、比較的単純に経度と緯度が求まります。まず、地球中心からGPS衛星までの距離をr、経度をλ、緯度をΦとすると、前回求めた、座標は以下のように表すことができます。

X_ECEF = r*cosΦ*cosλ
Y_ECEF = r*cosΦ*sinλ
Z_ECEF = r*sinΦ

これはいたって単純です。分からない人は、以前、お話したように、一度、地球とGPS衛星の絵を描いてみて、経度、緯度をベースに

さ～いん、こさ～いん

と指でなぞってみるとよくわかると思います。

さて、ここで求めたい値はΦとλです。まずは、経度λを求めて見ましょう。

まず、以下の式をみてください。


このように、まずはY_ECEF/X_ECEFが、tanλになることを利用し数式1の②でその逆数をとることで、経度λが求まります。注意せねばならないのは、アークタンジェントの結果は -π/2からπ/2の範囲でしか答えがえられませんので、求まった答えの象限に応じ加工してやる必要があるのですが、javascript等の言語ではatan2(Y,X)という感じで-πからπの範囲に変換してくれる関数がありますので、これを利用すると便利です。

次に、緯度Φを求めて見ましょう。地球を球と考えた場合、緯度も同様に求めることができます。


sqrt(X_ECEF²+Y_ECEF²)の部分は単純に3平方の定理で、Z_ECEFに対応するもう1辺を求めているだけです。あとは経度の時と同様に計算するだけです。いやはや、中学や高校のころは、三角関数がこんなところで利用できるとは思ってもみませんでした。

さあいかがでしたでしょうか？これで、経度緯度が求まったわけですが、実際には地球は楕円です。ですから、今回求めた方法よりも、計算は複雑になります。しかしながら、考え方や式の成分は今回と同じですから、とっつきやすいとは思います。
とはいえ、近似式を利用するなど、式自体を理解するには、軌道計算同様、時間を要します。

いずれにしましても、基本をまずはしっかり押さえておくことが大事です。

また、機会があれば、WGS84 準拠楕円体を考慮した、経度緯度の求め方も公開します。そのまえに、エフェメリスデータによる、GPSの位置計算もしたいところです。

[0回]

今回は、「WGS84測地系モドキによるGPS衛星の経度緯度の求め方其の1」でお話したWGS84測地系の座標系に変換してみます。

実は、ここの概念が、本件では一番分かりにくい話になります。

まず「アルマナックデータによるGPS衛星の位置計算其の3」で求めたGPS位置(X,Y,Z)の座標系はそもそも何だったのかからお話しますと、地球中心慣性ECI(Earth-Centered Inertial)座標系というものであり、座標原点が地球の質量中心にあり、ニュートン力学の原理基盤となっている座標系です。まずこの時点で挫折しそうですねえ。

次に、ちょっとすっとばして、上記の座標系を今回の主題であるWGS84測地系の座標系に変換するわけですが、この座標系は地球中心地球固定ECEF(Earth-Centered Earth-Fixed)座標系と呼ばれ、GPS 受信機などの位置を記述するために、地球の自転と共に回転する座標系です。

どうでしょうか？なんとなくわかってきましたか？まず、単純に言いますと、今回おこなう作業は地球の自転を考慮し座標系を回転させてやるのです。

ええ！じゃあ今までは、そうじゃなかったの？と思う人もいるかもしれませんが、そうじゃなかったのです。物事には順序というものがあります。まあ、がんばってついてきてください。


まずは図1をみてください。以前、x軸方向の座標を決定する際、春分点方向という言葉を使いました。実は、その時はさほど重要ではなかったのでさらっと流したのですが、実は、今回の座標変換の種を明かせば、この春分点方向から、地球の自転軸で回転させ、地球のグリニッジ子午線に合わせるという行為なのです。

グリニッジ子午線は経度0ですので、ここにあわせるのは緯度、経度を求める上では明快だと思います。しかしながら、そもそも春分点とは何か？という話なのですが、これは

「黄道と天の赤道との2つの交点（分点）のうち、黄道が南から北へ交わる方の点（昇交点）のこと。この点が赤経0 度かつ黄経0度であり、この点を太陽が通過する瞬間が春分となる。」

というものです。ますます意味不明ですよね。とりあえず、天文学的に、星の位置を示すための基準点とでも思ってください。

ここで、図1のθはグリニッジ恒星時とよばれている角度になるのですが、これを求めるためGPS衛星の場合は、航法メッセージにも定義されている、地球の平均自転速度ω_E=7,292,115.1467E-11を利用します。

実際にθを求めるには、GPS週の最初を、春分点方向 = グリニッジ子午線とし、そこから、現在時刻までの自転による角度のずれをω_Eを利用し求めるだけです。計算的には以下の通り。

θ = ω_{E * GPS週の最初から現在までの秒数}

すでにおわかりでしょうが、GPS週はアルマナックデータより取得可能です。尚、グリニッジ恒星時も奥が深いのでここでは深く突っ込みません。

ここでθがもとまりましたので、すでに求めている。GPS衛星位置(X,Y,Z)を回転行列を用いて、θ回転させてみましょう。


これはいたって明快ですよね。正直、すごく単純に記載しましたが、天文学も奥が深く、簡単に語れるものではありませんが、計算はこんな感じです。

今日はここまで。

[0回]