前回「アルマナックデータによるGPS衛星の位置計算其の1」で、アルマナックデータからGPS衛星の位置を求めるための材料をそろえました。さて、今回は離心近点角をケプラーの方程式を利用して求めます。

すでに前回、平均近点角Mを求めていますから、簡単に求まるようにも思えますが、残念ながら、すんなりとは求まりません。理由としては、非線形だからです。したがって、ニュートン法を用いて、近似解を求めます。尚、ここではニュートン法の説明はしませんのでがんばってググってください。

さてまず、ケプラーの方程式をニュートン法にて離心近点角Eを解くため以下のように設定します。


次に、以下の漸化式に数式1をあてはめます。


数式1を数式2に当てはめると次式が得られます。


ここで、初期E₀ = M = 73.40211084258863として、数式3を繰り返し演算していけば、値が収束します。永久に計算をし続けるわけにもいけませんので、前後のEの値を比較して、十分差が小さければ、そのEを利用します。

尚、ニュートン法では、解になるべく近い初期値を設定しないと、発散したり、なかなか、解が求まらないという自体になりかねませんが、Mを初期値にしておけば、この場合問題ありません。

さて、実際にEを求めるとE = 73.40125940150486 となります。

ニュートン法は軌道計算にくらべたら、まだ理解しやすい方なので、ググって是非理解してください。解析には欠かせないアイテムの一つです。

では次回で最終章といたしましょう。

[0回]

前回「航法メッセージによる GPS衛星の位置の求め方 其の3」で航法メッセージからGPS衛星の位置を求める概念をお話いたしました。さて、ここでは、航法メッセージの中に含まれるアルマナックデータを利用して、実際にGPS衛星の位置を計算したいと思います。

まずは、なるべく新しいアルマナックデータを取ってくる必要があるのですが、これは、米国政府が運用している、ナビゲーションセンターのサイトから入手可能です。

また、今回はまだ運用が開始していませんが、折角ですので、5月28日に打ち上げられた、ブロックIIFシリーズのGPS衛星の位置を計算してみます。

取得したアルマナックデータの内容は以下の通りになります。ついでに、変数も定義しておきます。

変数	意味	値
ID	衛星のPRN 番号	25
HEALTH	衛星の健康状態（000 の場合は利用可能）	063
e	軌道離心率	0.9350776672e-3
ta	現GPS週におけるalmanac軌道要素の基準時刻（s）	61440
i	軌道傾斜角（rad）	0.9600639343
Ω ⋅ dot	昇交点経度の変化率（rad/s）	-0.7425114745e-8
√a	軌道長半径の平方根	5180.759277
ℓ0	基準時刻における昇交点経度（rad）	-0.3115975261e+1
ω	近地点引数（rad）	2.185029387
M0	基準時刻taにおける平均近点離角（rad）	0.1111893773e+1
a0	時計の位相バイアス（s）	0.9536743164e-6
a1	時計の周波数バイアス（s/s）	0.0000000000e+0
WEEK	GPS 週	564

数字をみるとなにやら気分が悪くなりそうですね。前回の概念があれば、なんとなくは分かると思います。

計算に入る前に捕捉すると、IDはGPS衛星の番号です。そして、HEALTHは健康状態で、000以外の数字が入っていれば、利用不可能です。今回はまだ運用が始まってませんので、000以外の値が入っております。また、a0,a1は位置計算では利用しません。これは、受信機の内部時計の補正等に利用します。

さて、ここから本格的に計算を始めますが、まずは下ごしらえをします。計算に使う材料を計算させるとでも思ってください。

まず、平均運動N0を求めます。単位時間あたりにどれだけ衛星が地球を公転するかがこれを求めればわかります。これを求めるには、ケプラーの第3法則とニュートンの公式から、


と表すことができますので、GMは前回お話したとおり、既知でGM = 398600500000000、aはアルマナックデータの二乗で a = 26840266.686221562
になり、N0 = 0.00014357825789904194となります。

次が重要ですが、アルマナックデータが生成された、基準時刻と、GPS衛星の位置を求める時刻との差を求めます。これがないとはじまりませんよね。
まず、アルマナックデータが生成された基準時刻を求めますが、これはアルマナックデータのtaとWEEKを用います。

WEEKは、1980年1月6日0時0分0秒(UTC)を第0週として、取得したアルマナックデータを生成されたのが第何週目にあたるかを表しており、taはその週を基準とし、何秒後にアルマナックデータが生成されたかを表しております。

1980年ってのはGPSの運用の始まりと推測はつくと思います。ちなみにこの時点では4基のGPS衛星しか存在してなかったようです。

さて、ここまで言えば、時間差は容易にもとまるといいたいところですが、WEEKをよく見てください。どうみても、1980年からのWEEKにしては少なすぎます。実は、古い、GPS衛星のWEEK用のカウンターは10ビットつまり、1023までしかカウントできなかったため、1999年8月22日に1度リセットがかかっているためです。これは俗にゆう2000年問題の一つだったようです。ここらへんの対策はすでに、カーナビ等には備わっていたのか、問題はさほどでなかったようです。新しいGPS衛星では問題ないのでしょうが、今更、これを変更させるわけにもいかないので、今後どのようになるのか見ものです。

さて、というわけで、実際のWEEKは1024 + 564になります。

よって、1980年1月6日0時0分0秒からアルマナックデータ更新までの通算時刻(基準時刻(s))は

(7.0 * 24.0 * 60.0 * 60.0) * (1024 + 564) + 61440秒となります。

したがって、1980年1月6日0時0分0秒から現在の時刻(衛星の位置を求める時刻)の通算秒からこれを引けば、アルマナック生成の基準時刻からどれだけの秒数がたったのか分かります。ここで、(1980年1月6日0:0:0から現時点までの通算秒) - ((7.0 * 24.0 * 60.0 * 60.0) * (1024 + 564) + 61440) = 503490 = tkとします。

だいぶ材料がそろってきました。

さて、ここで、ようやく、現在時刻における、平均近点角を求めることができます。平均近点角Mは、以下の通り求めることができます。

M = M0 + tk * N0

M = 1.111893773 + 503490* 0.00014357825789904194 = 73.40211084258863
ここでM0はアルマナックデータより得られる。アルマナックデータ生成時の平均近点角ですので、そこを基準にどれだけ回転したかを求めるだけですので、いたってシンプルですよね。

さて、この時点で昇交点経度も求めておきます。残念ながら、昇交点経度も時間とともに変化しますのでこれも必須です。考え方は平均近点角の求め方と同様で、現在時刻による昇交点経度をΩKとすれば、

ΩK = ℓ0 + tk * Ω・dot

で求めることができます。したがって、Ωk =  -3.115975261 +   503490 *  -0.7425114745e-8 = -3.11971373202296
となります。

これで、仕込みは完了です。次回は、ケプラーの方程式等を用いて、本格的な計算をします。

[4回]

「航法メッセージによる GPS衛星の位置の求め方 其2」で2次元空間上でのGPS衛星の位置の求め方をお話いたしました。さてようやくですが、3次元空間上でのGPS衛星の位置を求めてみましょう。

まず、地球上のGPS受信機により、衛星位置は利用されますので、基準座標をGPS衛星の楕円軌道上の中心から地球の中心に移動させ、以下の図のようなイメージで座標空間を設定します。


図1は地球を赤道面で真っ二つに割ったイメージで、赤道面中心から、春分点方向にx軸、赤道面中心から北極点方向にz軸をとります。この時、必然的にy軸方向は決定します。

さて、ここで、前回お話した、楕円上の衛星位置(x',y')を求める式を以下の通り、地球中心とした形に変更します。


これに関しては問題ないと思いますが、x'とy'とダッシュをつけたのは、3次元での軌道でのx,yとは異なるためです。
まずは次の図を見てください。


なんかややこしい図がでてきましたが、まず頭をやわらかくして、聞いてください。図2に軌道面とあります。これは、GPS軌道の平面を表します。そして、この軌道面ですが、3次元上では、図のように、赤道面からみて傾いています。

ここまで、くるとなんとなく分かりますが、2次元上のGPS衛星の位置(x',y')について軌道面を斜めに向けた時の位置にしてやればいいわけです。

さて、3次元に座標を変換するための軌道パラメータが与えられていなければなんにもすることができませんので、ここで、簡単に航法メッセージより得られる軌道面に関する軌道パラメータを、紹介すると、

昇交点赤経Ω、軌道傾斜角i、近地点引数ωがあります。これはすべて角度をあらわし、図2に記述している通りで、これが決定すると軌道面の傾きがわかるわけです。とはいえ、なんかぴんとこないので、次のように考えると分かりやすいです。

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

こうすると感覚的にわかりますよね。よって回転行列を用いて、以下のようにGPS衛星の3次元座標上の位置(x,y,z)を求めることができます。


2次元面ではz軸は存在しませんので0になります。尚、ここでは回転行列の説明はしません。ただ、回転行列の概念は3Dをやる上では結構あそべるネタです。

一応、基本的な求め方はこれで終了なのですが、正確なGPS衛星の位置を計算するための軌道パラメータ要素は実際にはまだあり、これだけの計算では不十分と言えます。この話はまたの機会とします。

[0回]