IGSの放送暦を使った測位衛星の位置計算 |
IGSの放送暦を使った測位衛星の位置計算 |
IGSの放送暦を利用して測位衛星位置の計算をやってみました。 プログラムのコア部分には電子航法研究所（ENRI）の福島さんと坂井さんが入門用に公開されておられる GPS測位計算プログラム を利用させて頂いています。 このページは「GPS測位計算プログラム入門をやってみた」の続編と云う位置付けです。 プログラムの開発にはMicrosoft社が無料でリリースしているVisual C++ 2010 Express Editionを利用しています。
IGSから得られる放送暦とそのフォーマット 放送暦の利用アルゴリズム プログラムのダウンロード 精密軌道暦と比較
放送暦はいくつかのFTPサーバよりダウンロードが可能です。 FTPサーバについてはIGSのサイトに案内があります。 リンク先のページの内、Broadcast欄の一番右にある、CDDIS(US-MD)，SOPAC(US-CA)，IGN(FR)をクリックするとアドレスの案内があります。 その内CDDISのサイトのアドレスは、「ftp://cddis.gsfc.nasa.gov/gps/data/daily/ｙｙｙｙ/ｄｄｄ/ｙｙn/」となっています。 ただし、ｙｙｙｙは年，ｄｄｄは1月１日からの通算日，ｙｙは西暦の下2桁を表しています。 例えば、「ftp://cddis.gsfc.nasa.gov/gps/data/daily/2010/010/10n/」です。 ダウンロードはブラウザからマウス使ってポチポチと手動でも可能ですが、専用のスクリプトが有ると便利です。
ファイル名は"brdcｄｄｄ0.ｙｙn"となっています。 ｙｙは西暦の下2けたを表し、ｎはNavigationのｎを表しています。 ただし、ファイルは圧縮されていますので、解凍（展開とも云う）が必要です。 お使いのマシンがWindowsであれば7-zipをお勧めします。
IGS(International GNSS Service?)が提供している放送暦を利用して衛星の位置を計算するために、放送暦のフォーマットを押さえておく必要があります。 IGSの提供する放送暦は図1.2.1に示す様な、RINEXというフォーマットになっています。 このRINEXフォーマットのデータは数年毎に細かい改訂が行われており、データの表記の仕方がバージョン毎に少しずつ異なります。 例えば、ここに示したものより古いフォーマットでは整数部分の"0"が省略されるなどです。
RINEXフォーマットの詳細についてはIGSのサイトを参照してください。（http://igscb.jpl.nasa.gov/components/formats.html）
2 NAVIGATION DATA RINEX VERSION / TYPE
CCRINEXN V1.6.0 UX CDDIS 20-JUN-10 02:51 PGM / RUN BY / DATE
IGS BROADCAST EPHEMERIS FILE COMMENT
0.5588E-08 0.1490E-07 -0.5960E-07 -0.1192E-06 ION ALPHA
0.8397E+05 0.9830E+05 -0.6554E+05 -0.5243E+06 ION BETA
-0.186264514923E-08-0.710542735760E-14 61440 1589 DELTA-UTC: A0,A1,T,W
15 LEAP SECONDS
END OF HEADER
1 10 6 19 0 0 0.0-0.132149551064E-03-0.397903932026E-11 0.000000000000E+00
0.490000000000E+02-0.360312500000E+02 0.483591572101E-08 0.288395691448E+01
-0.193342566490E-05 0.481024384499E-02 0.534579157829E-05 0.515479243088E+04
0.518400000000E+06-0.298023223877E-07-0.310813213882E+01-0.119209289551E-06
0.965296028888E+00 0.279937500000E+03 0.877006962495E+00-0.818534095231E-08
0.367872466222E-10 0.100000000000E+01 0.158800000000E+04 0.000000000000E+00
0.200000000000E+01 0.630000000000E+02-0.190921127796E-07 0.490000000000E+02
0.512880000000E+06 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00
2 10 6 19 0 0 0.0 0.265788286924E-03 0.329691829393E-11 0.000000000000E+00
0.550000000000E+02 0.570625000000E+02 0.500878005738E-08 0.122865347467E+01
0.315718352795E-05 0.957962975372E-02 0.925920903683E-05 0.515359644890E+04
0.518400000000E+06 0.521540641785E-07-0.102519684372E+01 0.162050127983E-06
0.939550794850E+00 0.188562500000E+03 0.308919237482E+01-0.826641599616E-08
-0.478591367925E-10 0.100000000000E+01 0.158800000000E+04 0.000000000000E+00
0.200000000000E+01 0.000000000000E+00-0.172294676304E-07 0.550000000000E+02
0.511218000000E+06 0.400000000000E+01 0.000000000000E+00 0.000000000000E+00
to be continued
図1.2.1 RINEXフォーマットによる放送暦の例
プログラムは放送暦の記録ファイルを読み込んだ後、RAM内に放送暦の情報を数値化して保存します。 その際、プログラム作成上の都合により読み込まれたデータはファイル先頭から順番に格納されます。 このデータを利用する上では、エフェメリスの放送時刻（TOM: Transmittion time of Message）に注意を払う必要があります。
IGSから得られる放送暦は実際のデータを見て頂くと分かるように、エポック順に並んでいます。 ここで云うエポックとは、記述された軌道情報の基準となる時刻の事です。 実際、図1.3.1に示す様にTOMに注目するとTOMが前後してしまっている事が分かります。 放送暦が放送時刻順に並んでいないため、RAMに蓄えた放送暦を格納した順で使うとTOMの古いものを使ってしまう可能性があります。
「GPS測位計算プログラム入門」がデフォルトで対応している処理は、ある受信機独自のフォーマットによる生のエフェメリスデータ（バイナリの16進表現）に対するものです。 従ってIGSのエフェメリスを読み込んで衛星位置算出に使えるようにする必要があります。 このプログラムは衛星位置の計算プログラムにENRI製計算エンジンを利用した上で、IGSより取得したGPSの放送暦（エフェメリス）を使った衛星位置計算を行います。 対応確認済みのRINEXバージョンは2002年以降のものですが、文字処理部分を工夫していますのでそれより以前のものでも構わないと思います。
このプログラムの第一の目的はエフェメリスの精度評価です。 従ってプログラムの計算の仕方や引数の与え方などはそれに最適化されています。 本プログラムを使用して求まるものは、あるエポック（任意）における衛星の重心位置座標（ECEF系）です。 IGSがリリースしている精密軌道暦（重心位置）と比較が可能です。
オプション的ものですが、非常に長い期間の計算を一つのエフェメリスで行いたい場合の改造点について説明します。 ENRIの衛星位置計算プログラムでは、引数で渡されたエポックとエフェメリス情報のエポックとの差を計算における変数としています。 デフォルトのプログラムでは、この差の絶対値が一週間の604,800[s]の半分である302,400[s]を超えると一週間前後にずらします。 従って、同一のエフェメリスを使って1週間を超える計算を行う場合はソースコードの改造が必要です。
[2010/8/24 追記] と云うか、それ以前にエポック2.5～が3日程度経過すると計算結果が不正になる問題が有ることが分かりました。 どうすればいいのか、今の段階では分かりません。位置計算のアルゴリズムは本家の「GPS測位計算プログラム入門」に譲るとして、ここではIGSの放送暦の利用に関するアルゴリズムについて説明します。 このプログラムの基本的な流れを図2.3.1と図2.3.2に示します。 先ず、引数として計算期間を受け取ります。 実行時に引数を受け取れなかった場合、scanf関数を使って計算期間の入力を求められます。 その後に計算ルーチンに入ります。 細かい点は以下のフローチャートを見て頂くとして、使用する上で注意の必要な事項について説明します。
計算に使う放送暦情報は、計算しようとするエポックより以前で放送されたデータの内最近のものを使います。
前述したTOMの不連続性に対処するため、デフォルト状態での計算ルーチンでは計算しようとしているエポックよりも約6時間分程度までの未来の情報が必要となります。 また、初回の計算のためには前日分のデータも必要となります。 従って、ある期間の計算を行おうとした場合、前後丸1日分の放送暦の用意も必要となります。
2010/11のGPS/GNSSシンポジウム2010でENRIの坂井さんより実行形式ファイルを公開しても良いと快諾頂きました。 この場でお礼を申し上げます。 と云う訳で早速公開致します。 折角ですので坂井さんの本を宣伝したいと思います。 このプログラムを開発するにあたって度々参考にしている書籍で、坂井 丈泰 (著) 「GPSのための実用プログラミング」です。 Amazonで\3,465です。 ぜひ中古ではなく新品を購入しましょう。
| リリース日 | プログラムパッケージ | 更新内容 |
| 2014/2/28 |
放送暦を基に衛星軌道を計算 with Python 使用言語：Python |
エフェメリスの選択アルゴリズムが全く異なるのですが、参考までにPythonで実装した衛星座標の計算スクリプトへのリンクを張っておきます。
以後、GNSSに関するプログラムのほとんどはPythonで作る予定です。
開発環境：Windows 7 (64/32 bit) |
| 2013/10/20 |
gnss_sat_position_VS2012.7z 使用言語：C++ |
バグの修正，準天頂衛星（QZS）への対応，出力ファイルフォーマットの変更，IGSの放送暦ファイルのフォーマットエラーへの対応（IGSの放送暦生成スクリプトにはバグがあると思う）を実施しました。
圧縮ファイルは、ソースコード，ビルド済みの実行形式ファイル，実行ファイルを簡単に利用するためのPowerShellスクリプト，データ処理用のPythonとRubyのスクリプト，説明書を含んでいます。
今回同梱している実行用のスクリプトはbatファイルではなく、PowerShellに変更しています。
これでパラメータの変更を行いやすくなっています。
また、出力されたファイルは、Rubyスクリプトを用いて処理することで、u-centerで衛星軌道を確認できるようになりました。
なお、今回のプログラムから、一部にC++の機能を利用しています。
開発環境：Windows 7 (64 bit), Visual Studio 2012 pro 開発環境が更新されていますのでご注意ください。  u-centerで描画した、指定座標から見た測位衛星の軌跡 |
| 2011/2/11 |
sat_position_calculation_for_release_v2.7z 使用言語：C++ |
デフォルトの状態で、直ぐに計算が行えるようにバッチファイルの定数を調整しました。
また、既に計算した結果も含めていますので、実行前にどんなファイルが作成されるのか確認することも可能です。
開発環境：Windows7 (32bit), Visual C++ 2010 Express Edition |
| 2010/11/9 |
SatPosCalcForRrelease.zip 使用言語：C++ |
初公開。 2010年の放送暦の一部を添付しています。 使用の際には放送暦は解凍の上で 実行形式ファイルと同じディレクトリに置いて下さい。 この圧縮ファイルは、計算に必要なコンパイル済みの実行形式ファイルと私＆後輩で作ったソースコードを同梱しています。 この後の章での説明に使用しているプログラムより更新した物ですので、ファイルの出力状況が若干異なります。 [2010/11/19追記]リンクが張れていなかったのを修正しました。 [2011/1/10 追記]ファイル容量制限に引っ掛かって定期的に削除されていたので同梱する放送暦を減らしました。 開発環境：Windows7 (32bit), Visual C++ 2010 Express Edition |
計算の実行は非常に簡単です。 先ず図2.5.1に有る様にバッチファイルをダブルクリックするだけです。 バッチファイルを実行後、図2.5.2の様に計算エポックが記入されたcsvファイルが出力されます。 その後、図2.5.3のに示す様にコマンドラインで衛星位置が計算されます。 最後には図2.5.4の様にlogフォルダ内に結果が出力されます。
計算パラメータは、計算開始日・計算終了日・インターバル（計算の間隔）・オフセット時間の4つです。 これらのパラメータを変更するには、図2.5.5に示す様にバッチファイルをテキストエディタで編集すればOKです。 図2.5.5の例では、2010/7/1から2010/7/4までの間で3600sec毎に計算することになります。
[注意1]：本プログラムは計算期間の前後1日分を含む放送暦が必要です。
もし、ファイルが存在しない場合でもエラーは出力されなかったと思います。
[注意2]：バッチファイルで指定できるインターバルは24時間以下である必要が有ります。
2010/11/9付けでリリースしたプログラム（2.4節参照）を利用して、計算を行います。 ここでは、IGSのエフェメリスと比較がやり易いように、900秒毎にエポックを切り替えながら計算を行います。 図3.1.1に示す様に、計算期間を2010/9/8から2010/9/11とし、インターバルを900秒にセットしました。
プログラムを実行すると、図3.2の様に計算結果がcsvファイル形式で出力されます。 「prn0_log.csv」には、全衛星のデータがエポック順に並んで出力されています。 その他のファイルは個々の衛星の位置情報です。
標準で出てくる衛星座標はECEF座標系ですのでx,y,z表示です。 これを手っ取り早く表示して見せるためには、Rコマンダーを使うのが近道の一つです。
次に、地上との関係を知りたいと思ってGoogle Earthで表示させたのが図3.2.2です。
出力されるx,y,z座標を緯度・経度・高度に直して、kmlフォーマットに成形すると容易に表示が出来ます。
地球をぐるぐる回して動かせるので面白いです。
ここで公開しているプログラムは衛星の座標をx,y,zで出力しますので、これをBLH（緯度・経度・標高）へ変換してさらにKMLフォーマットにしなくてはなりません。
これは結構面倒ですのでスクリプトを用意しました。
下記ページよりダウンロードして下さい。
ただし、BLHへの変換には近似式を使って、標高は楕円体高で代用しています。
[x,y,zをBLHへ変換するスクリプト]
[BLHをKMLへ変換するスクリプト]
3章での計算で得られた結果を、IGSが公開している精密軌道暦と比較します。 精密軌道暦はIGSの「http://igscb.jpl.nasa.gov/components/prods_cb.html」よりダウンロードすることが出来ます。 GPS週番号をクリックして、その中から「igswwwwd.sp3.Z」をダウンロードして下さい。 ただし、wwwwはGPS週番号であり、dは曜日番号（日曜日が0）です。 ダウンロード後は放送暦をダウンロードした時と同様に解凍と拡張子から「.Z」を削除する作業を行って下さい。 精密軌道暦は計算結果の格納されたフォルダの中に、ディレクトリを作ってそこに格納しておきます。
図4.1.2に示す様に、「estimate_all.rb」と云うRubyのスクリプトを実行することで精度評価を行います。
このスクリプトについては以下よりダウンロードして下さい。
[精度評価スクリプトのダウンロード]
このスクリプトを実行すると図4.1.3の様に、「prn*_log_error.csv」と「eph_error_for_R.csv」と云うファイルが生成されます。 「prn*_log_error.csv」は個別の処理結果が格納されており、「eph_error_for_R.csv」はR（アール）と云う統計処理ソフトで読み込みやすい様に成形されたファイルになっています。
精密暦と比較した結果を図4.2.1～図4.2.4に示します。 横軸が時間で縦軸が各軸の誤差を表します。 図4.2.1には2010/7/7-2010/7/10におけるPRN2の衛星の軌道予測誤差をx,y,z成分で表示しています。 計算に使用される放送暦が切り替わるたびに計算値に差が出ているようですが・・・これで良いんですよね?
次に、全衛星の結果を示します。 PRN番号と衛星のシリーズを比較すると、新型で有ればある程精度が良いようです。 また、同じブロックⅡRM衛星でも、打ち上げ時期が最近の物ほど精度が良いように見えます。 これは搭載周波数源の原子時計の性能向上によるものでしょうか…。
細かく見ると、殆どの衛星の予測誤差には減少又は増加の傾向が見られます。
3日間の結果を見る限り、周期の長い変動が見られそうでしたので思いっきり長い期間の計算を行ってみました。 ここではPRN2とPRN8について図3.4と図3.5に示します。 この場合、たまたま同じ様な変動を示しています。 変動の周期を計算してみると月の満ち欠けの周期に近い物となっていました。 月と何か関係有るのでしょうか?
こちらはおまけで同じ放送暦を使い続けた場合にどの様に軌道の予測精度が劣化していくのかを計算した結果です。 基本的には、基準エポックから2～4時間経過すると誤差が急速に増大する様です。 また、最近の衛星であるほど誤差の大きさが小さい様な気がしました。
