Pythonを使ってGPS,QZSS現在位置表示（軌道要素編) - やべべのべ(雑用系リーマンの独り言)

個人的にはニッケイ、ダウのソウバが気になってしょうがないのですが、
3か月ぶりにGPS関連です。前回の続きです。↓前回

tamhatam.hatenablog.com

こちらの記事はニーズがないので誰も見てないと思いますが、

備忘録もかねて記載していきたいと思います。

Pythonでの初アプリとなるGPS/QZSSビューワですが、
今回は基幹部分の軌道計算について記載したいと思います。
（書いてて１回では終わらなさそうです。）

まず、衛星の位置を知るためには、衛星の軌道を知らなけらばなりません。
超基本的には、”衛星は地球の周りをぐるぐる回っているだけ”ですので、
軌道円と時刻がわかれば、位置座標が決定されるわけです。
回る速さについては、衛星の高度できまります。
GPS/QZSSはとても単純でGPSは約12時間、QZSSは約24時間で地球の周りを
１周しています。

GPS/QZSS共に衛星軌道は、
だいたいの軌道を表現した「アルマナック」と
精密な軌道を表現した「エフェメリス」に区別されます。

エフェメリスは、受信機内部で使用されており、自己位置算出に用いるものです。
対してアルマナックは衛星のだいたいの位置を予測する目的で使用するものです。
受信機では、アルマナックを計算してどの衛星がどこにいるか？計算して、
捕捉したい衛星を狙う際に使用しているイメージです。

今回はアルマナックを用い、衛星の位置を算出してみます。
衛星の位置座標がわかれば、自分の位置(例えば東京)から衛星がどこに見えるか？
わかるわけです。東西南北、地平線からどのくらいのところに見えるのか？
といったところです。

アルマナックですが、QZSSのHPで毎日最新のものが手に入るため、
これをダウンロードし、読み込んで使用することとします。
こちらのフォーマットは以下となります。
まずファイル名ですがqg20200016.alm
です。これは"q"zssと"g"psのアルマナックで"2020"年の"016"日目に公開されたものですよという意味です。

んで、中身は以下です。

******** Week 41 almanac for PRN-01 ********
ID: 　　　　　　　　　　　01
Health: 　　　　　　　　　000
Eccentricity: 　　　　　　　9.284019470E-03
Time of Applicability(s): 　　61440.0000
Orbital Inclination(rad): 　　0.9785982499
Rate of Right Ascen(r/s): 　-7.657461821E-09
SQRT(A) (m 1/2): 　　　　　5153.597168
Right Ascen at Week(rad): 　-9.528350126E-01
Argument of Perigee(rad): 　0.756456635
Mean Anom(rad): 　　　　　1.845888966E+00
Af0(s): 　　　　　　　　　　-2.679824829E-04
Af1(s/s): 　　　　　　　　　-1.091393642E-11
week:　　　　　　　　　　　 41

１個づつ解説していきましょうか

ID:　運用している衛星の識別番号です。 GPSは1～32,QZSSは193～199(現時点)です。

Health:　衛星の状態を表します。000は正常です。000以外は基本使用できる状態ではない場合がほとんどです。例えば軌道修正等メンテや問題が発生して運用停止中などに000以外になります。よくみるのは"031"です。(詳しくはICDを参照。)

Eccentricity:　離心率です。

Time of Applicability(s):　このアルマナックの基準となる時刻です。UTC時刻で日曜日のAM0:00:00を0としています。この基準時刻から何秒後(前)なのか？=>軌道のどこの位置なのか?になります。

Orbital Inclination(rad):　軌道傾斜角です。楕円軌道が赤道面に対してどれだけ傾斜しているか?(rad)です。

Rate of Right Ascen(r/s):　昇交点赤経の時間変化率。赤道面と衛星軌道面が交差する点(ただし、南から北に上がるところ）春分点とのなす角度が昇交点赤経(Ω)でその時間変化率。参考URLを見てください。

SQRT(A) (m 1/2): 　軌道長半径の平方根です。

Right Ascen at Week(rad):　週の初め(週秒でいうと0)の昇交点赤経(昇交点経度といっているサイトもあった。）

Argument of Perigee(rad):　近地点離角(近地点引数とも)ω。衛星軌道面の近地点と赤道面とのなす角度。参考URL見て。

Mean Anom(rad): 　平均近点角です。

Af0(s): 　時計の位相バイアス（s）　→今回は使用しません

Af1(s/s): 　時計の周波数バイアス（s/s）　→今回は使用しません。

week: 　GPS週です。GPS週は10bitで表現されるため、0-1023までで表現されます。GPSは1980年1月6日午前0時を起点としており、1999年8月21～22日で1023→0に繰り上がり、2019年4月6日～7日で2回目の繰り上がりが起きました。(これをロールオーバといっています）長くなりましたが、2019年4月7日からの経過週数です。

長くなったので、今回はアルマナックの説明で終わっちゃったんですが、
このアルマナックという軌道要素を用いて衛星の位置を算出します。

流れ的には、

軌道長半径、離心率から2次元の衛星軌道（楕円)がが決まり、
平均近点角と基準時刻と位置を知りたい時間との差分から衛星軌道のどの位置にいるか？算出します。
衛星軌道とその軌道のどこに衛星がいるか？算出できたら、軌道傾斜角や昇交点なんちゃらを使って衛星軌道面と赤道面との関係から2次元を3次元にします。（楕円軌道をECEF座標に変換）

となります。

まあ文章ではよくわからんと思うので、以下サイトを参考としてください。
めちゃわかりやすいです。(プログラム作る際に参考にさせていただきました）

っていうか以下サイトめちゃわかりやすいので正直
全部書いてあるので自分のブログは読まなくてオケ。。。

アルマナックからGPS衛星の位置を求める方法（衛星位置算出全般）

アルマナックからGPS衛星の位置を求める方法 | Nyanchew's Digital Life - ja1

JAXA衛星セミナー第1回衛星の軌道と打上げ時刻（衛星軌道要素の説明）

http://www.satnavi.jaxa.jp/basic/seminar/pdf/seminar_090611.pdf

人工衛星の軌道（アルマナック、2行軌道要素形式、エフェメリス）（衛星位置算出全般）

人工衛星の軌道（アルマナック、2行軌道要素形式、エフェメリス）

4月7日（日本時間）に2度目の「GPS週数ロールオーバー」（ロールオーバのはなし）

qzss.go.jp

また、アルマナックとは別に衛星の軌道を表すものとしては、2行軌道要素形式(TLE)があります。というかこっちのが一般的な気がします。
アメリカの軍？のサイトに登録すれば、商用の衛星はほぼすべて情報が手に入るようです。