アルマナックからGPS衛星の位置を求める方法

GPS Test、GPS StatusなどAndroid端末のGPS信号受信状況を表示するアプリがいくつかあります。それらのアプリは現在のGPS衛星の位置を天空図に表示しています。そのためには現在位置からみた各衛星の方位角（Azimuth)と仰角（Elevation)が必要なのですが、これらのデータはアルマナックという衛星の軌道情報から計算で求めることができます。

自分でもやってみようと思って調べながらプログラムを書き、試行錯誤の結果、GPS衛星の位置を求めることができるようになったので、忘れないためのまとめとして書いておきます。

アルマナックというのはGPS衛星の軌道情報データで、衛星のおおよその位置を求めるために使います。今回はQZ-visionで配布しているものを使いました。Web APIでデータを取得することができます。

http://qz-vision.jaxa.jp/USE/ja/almanac

アルマナックの例を以下に示します。

******** Week 739 almanac for PRN-01 ********
ID:                         01
Health:                     003
Eccentricity:               2.453327179E-03
Time of Applicability(s):  147456.0000
Orbital Inclination(rad):   0.9602863543
Rate of Right Ascen(r/s):  -8.0117622935E-09
SQRT(A)  (m 1/2):           5153.626465
Right Ascen at Week(rad):   2.189035151E+00
Argument of Perigee(rad):   0.416642440
Mean Anom(rad):             2.076584921E+00
Af0(s):                     8.678436279E-05
Af1(s/s):                   3.637978807E-12
week:                        739

おおよその手順は以下のようになります。
1．現在時刻に衛星が楕円軌道のどの位置にあるかを地球を原点とする楕円平面の2次元座標系で求める
2．もとめた座標を地球の中心を原点とする3次元の座標系（地球中心地球固定座標系：ECEF）に変換する
3．さらに方位角と仰角を求めたい地球上の場所を原点とする地平直交座標（ENU座標）に変換する

衛星軌道上のどの位置にあるかを求める

GPS衛星は楕円軌道を描いて地球の周りを回っています。地球は楕円の中心ではなく焦点の位置にあります。

図1

楕円は以下の式で表されます。

…………(1)

ここで、a=長半径、b=短半径です。アルマナックには aの平方根(Square Root of Semi-Major Axis)が含まれています。地球は楕円の中心からa * eだけ離れています。e は離心率といいます。離心率eは以下の式で表されます。

…………(2)

まずはGPS衛星の楕円軌道上の位置を2次元で求めるわけですが、これに必要なパラメータや関係式がInterface Control Document という文書のTable 20-IV. Elements of Coordinate Systems にのっています。基本的にはここに書かれているとおりにやります。以下は現在時刻の衛星位置を求める手順です。

図1の平面上での地球を原点とした衛星のXY座標を求めるためには、図のV_k（真近点角）を求めればよいのですが、そのためにはE_ｋ（離心近点角）を求めることが必要です。つまり、衛星の楕円軌道上の地球を原点とした2次元座標は

…………(3)

となります。R_kは地球と衛星との間の長さで、以下で計算できます。

…………(4)

ここでV_kは離心近点角E_kを用いると以下のように書けます。

…………(5)

したがって、

…………(6)

となります。
E_kはどうやって求めるかというとここでケプラーの方程式が登場します。

…………(7)

M_kは平均近点角といってアルマナックに含まれる情報から計算できます。

…………(8)

M0はアルマナックに含まれる基準時刻における平均近点角(Mean Anomary)、T_kは現在時刻、N0は平均運動といって以下の式で求まります。GMは定数で 3.986005E+14です。

…………(9)

現在時刻Tkは、アルマナックデータの基準時刻(Time of Applicability)からの経過秒でなければならないので、GPS週の値とTime of Applicabilityを使って計算します。Time of Applicabilityは現在のGPS週の始まり（日曜日の午前0時）からの経過秒なので、以下の式で求めることができます。

…………(10)

T_gps は現在時刻TkをGPSエポック(1980年1月6日0時0分0秒）からの経過秒で表した値です。また、T_oeはアルマナックの基準時刻(Time of Applicability)です。week はGPS週番号です。GPS週番号に1024足しているのは、1999年8月にGPS週が1023を越えて一度リセットされているためです。

以上から(8)の式を解いてE_kを求めます。(8)はニュートン法を使って解きます。

…………(11)

とおいて、

…………(12)

にあてはめると、

…………(13)

となるので、初期値 E₀ をM_kとして上式を繰り返し計算しある程度Eの値が収束したところで、その値をE_kの値とします。

(4)からV_kを求めて(3)から座標を求めれば楕円軌道平面上の2次元座標系における座標値X’,Y’を求めることができます。

楕円軌道座標をECEF座標に変換する

次に求めた楕円軌道面上の2次元座標を地球中心の3次元座標（ECEF:Earth-Centered,Earth-Fixed)に変換します。

ECEF座標系は地球の中心を原点とし、北極方向をZ軸、グリニッジ子午線方向をX軸とする3次元の右手座標系です。衛星の楕円軌道は天球に対して固定であるのに対して、ECEF座標系は地球の自転といっしょに回転しています。

さて、楕円軌道面とECEFの赤道面（XY平面）との関係から座標系の回転を行ってECEF座標を求めます。

ここのところがわかりにくいので、別記事を書きました。

衛星の楕円軌道面と地球の赤道面との関係

まず、V_kにアルマナックの近地点引数(Argument Of Perigee)ωの値を加えた角度P_kを求めます。

…………(14)

(3) で V_kの代わりにP_kを使ったX’,Y’を求めます。

次に、軌道傾斜角 i (Inclination)だけX軸まわりに回転させ、さらにΩ_kだけZ軸まわりに回転させればOKです。i はアルマナックに含まれていますが、Ω_kは以下の式で計算します。

…………(15)

Ω_eは地球の自転速度で7.2921151467E-05です。Ω_dは基準時刻における昇交点黄経(Longitude of Ascending Node of Orbit Plane at Weekly Epoch)です。 Ω₀は基準時刻における昇交点経度変化率(Rate of Right Ascention)です。T_oeは基準時刻（Time of Applicability)です。

以上から、ECEF座標系のX,Y,Zは以下の回転行列で求めることができます。

…………(16)

これを展開すると以下のようになります。

…………(17)

ECEF座標をENU座標に変換する

最後にENU座標に変換します。ENU(East, North, Up)座標は地表面のある点を原点として天頂方向をZ軸、東方向をX軸、北方向をY軸とした座標系です。

ECEFからENUへの変換は座標系の回転と原点の移動でできます。ENU座標の原点にする地表面の位置の経度をλ度、緯度をφ度とすると、

・Z軸中心にλ度回転する（Z軸の正方向に向かって右ねじの方向に回す。グリニッジ子午線のところにあった基準子午線が原点の子午線のところになる）
・Y軸中心に(90度ｰφ度)回転する（Z軸が原点の天頂方向を向くようになる）
この時点でY軸が東、X軸が南を向いているので、
・Z軸中心に90度回転する
・原点を移動させる

式で書くと以下のようになります。

…………(18)