GPS測位計算プログラム入門をやってみた～VC++への移植～ |
電子航法研究所（ENRI）の福島さんと坂井さんが入門用に公開されておられる GPS測位計算プログラム を使ってみました。
Linux版（C言語）をコンパイル Windows版（Borland C++）をMicrosoft VC++でコンパイル VC++への移植上の問題 精密軌道暦と比較
「測位計算プログラム入門」とは、電子航法研究所（ENRI）の福島さんと坂井さんによって公開されているGPSの測位演算原理の学習プログラムです。 プログラムは3ステップに分かれており、それぞれ以下の様なテーマが設定されています。
その3では衛星位置が既知としてして計算が行われています。
この衛星位置についてはその2で準備を行います。
という訳で、このページではそんなサンプルプログラムを通して色々試した結果を備忘録的に記録しています。
まずはプログラムとエフェメリスデータファイルをダウンロードします。 今回はLinux系OSであるUbuntuマシンで実行したのでLinux版をダウンロードしました。 ちなみに、Windows OSにVirtualBoxをインストールして使えるUbuntuのVirtual Box用仮想マシンでも正常なコンパイルと実行を確認しています。
ダウンロードが出来たら圧縮ファイルを展開します（Windowsで言う解凍）。
後はコンパイルがすぐにでもできるのですが、一応必要なパッケージがあるかチェックをします。 パッケージマネージャを開いて確認しました。
build-essential がインストールされていればOKです。 インストールされていない場合は項目名をダブルクリック後に適用をワンクリックでインストールが開始されます。
では準備が出来ましたので実際にコンパイルをやってみます。 コンパイルするディレクトリ（フォルダ）の中で右クリックして端末を開きます。
"make"と入力してエンターキーを押すとコンパイルが実行されます。
次に、2nd（その2）の方もコンパイルを実行します。 下図の様にコンパイル時に警告が出ましたが、取り敢えず大丈夫です。
このまま実行させることも出来ますが、始めは1行しかデータの処理と表示を行いません。 そこで下図にあるように1行に制限するコード部分をコメントアウトしました。 この後はファイルを上書き保存して再びmakeを実行します。
↑before
after↓
予めrepa.datファイルを同じディレクトリに移した上で“./get_eph rapa.dat”と入力して実行すると下記の様に出力されました。 各衛星の軌道情報の一部の様です。
今回は実行環境にVC++ 2010を選択しました。 インストールと製品の登録が終わったら、プログラムをコンパイルする準備を行います。 ここで用いたのはBorland C++用のソースコードでしたが、たぶんLinux用のファイルでもコンパイルはできると思います（解凍は別にして）。
ビルドをやり易くするため、表示を上級者用に変更します。
適当な名前でプロジェクトを作ります。
開発プログラムのタイプはwin32コンソールを選びます。んで“OK”をクリックします。次の画面では“次へ”をクリックしてください。
“空のプロジェクト”にチェックを入れてから完了をクリックします。
ダウンロードしたプログラムを予めプロジェクトフォルダへ移しておきます。んでプロジェクトにプログラムソースを追加します。
プロジェクトにソースコードが追加されました。
ビルドの種類に“Release”を選択します。
エラーがなければ、この様に表示されます。
ビルドを実行するとプロジェクトフォルダ内に実行形式ファイルが生成されます。
出来たプログラムを実行するために、コマンドプロンプトのショートカットをコピペします。
そのままでは作業フォルダまでパスが通ってくれませんので、ショートカットのプロパティを開いて作業フォルダ欄を空にします。
これで準備完了です。コマンドプロンプトを立ち上げて実行させてみましょう。
プログラム名を入力してエンターキーを押すとこの様に実行されました。 座標変換が出来ています。 ちなみに、デバッグモードで動かすと分かるのですがPCで動かすために必要なプログラムが色々関連付けられています。 ですのでデバッグモードを使って変数の動きを見たい時は適当な位置にブレークポイントを置くことをお勧めします。
その1のソースコードは単体のソースコードでしたので殆ど問題が有りませんでしたが、その2では複数のソースコードをコンパイルする必要がありますので若干テクニックが必要です。 ソースファイルの登録方法には2つの方法があります。 一つ目は、図2.2-1に示す様にソースファイル欄に全て追加してしまう方法です。 これは“ソースファイル”の上で右クリックから「既存の項目の追加」を選ぶと出来ます。 2つ目は、図2.2-2に示す様にmain関数のあるeph2ecef.cに参照するソースファイルを全てインクルードしてしまう方法です。 ところで、この作業さえしてしまえば後はビルドするだけ…とはなりません。 2.3節で解説するエラーについて対処する必要があります。
 図2.2-1 ソースファイルとして登録した場合 |
 図2.2-2 インクルードする場合 |
その2をコンパイルしようとするとエラーや警告のオンパレードです。 この節ではそのエラー内容毎にどの様な方法で対処するか解説します。
警告の殆どは文字列処理関数やファイルの入出力関数のバッファオーバーランに関するセキュリティ関係です。 他には、使用していない変数が多数あったり、条件によっては返り値が不定となり得る関数があったりです。 また、他のソースコード中に宣言されたの関数を呼び出せなくて定義なしと出力されています。 さらに、グローバル変数のサイズが異なる変数宣言もエラーが出力されました。
strcpyやstrncpy，fopenはより安全なstrcpy_s等への変更が推奨されています。 これはバッファオーバーランを起こした場合、strcpyではそれを検出できないためです。 恐らく、やろうと思えばパソコンのRAMを引っかき回した揚句OSをダウンさせることもできるのではないでしょうか? strcpy_sはオーバーランの直接的な防止はしないらしいのですが、検出機能を持っているので危険を察知できるという訳です。 そこで、msdnを参考にして以下の宣言をeph2ecef.cのヘッダ部分に追加します。
将来的な拡張性を考えられているためか、そのままでは未使用の変数があります。 そのままにしておくとwarningが沢山出ますし、他のエラーが見えなくなるので全てコメントアウトすることで対処しました。
場合によって必要のない処置ですが、やっておいて損はありません。 各ソースコードでインクルードしているヘッダファイルを一度だけ確実にインクルードするための処置です。 C言語では定石の#ifdef, #define, #endifを使う方法です。 このコードを以下のコードの例のように各ヘッダに追加していきます。
tbufと言うグローバル変数が非常に重要な役割を果たすのですが、メインであるeph2ecef.cとread_data.cで宣言されているサイズが異なるためにコンパイラがエラーを出力します。 もし出力されていない場合でも実行させるとうまく動かないはずです。 という訳で、この変数はどちらか片方の宣言の仕方に統一することにします。 大きい方が安心なので今回はread_data.cの方を以下の様に修正しています。
extern char tbuf[FRM_MAX]; －＞ extern char tbuf[FRM_MAX+1];
eph2ecef.cは他のソースコードから様々な関数の呼び出しを行いますが、その関数が見つからないというエラーが出力されました。 この場合は、他のソースに関数を公開する意味でexternを使います。 関数を宣言しているファイルと、eph2ecef.cの先頭の各ヘッダ部分（上の方という意味）で外部関数宣言を行います。
その2のプログラムを実行すると衛星の位置座標が受信機内部の計算値と、観測された擬似距離を使った計算で求めた結果の比較が出力されます。
これを少し変更して、エフェメリスを使った衛星位置について精密軌道暦と比較してみましょう。
ここではIGSの観測データをこちらよりダウンロードしました。
その2ではサニャック効果を考慮した演算を行って、衛星の位置座標を求めています（ここがややこしい）。 衛星の位置は測距信号の発信された時刻で計算しますが、それを表す座標系は受信した時点のECEF座標系でなければなりません。 GPS受信機が記録した擬似距離と衛星座標の計算値は、あるエポック時刻において観測（計算）されたものです。 計算原理から言えばそのエポックで記録した擬似距離を光速で割って求まる時間分、測距信号は過去に発射されたことになります。 この発信時刻の推定精度を上げるために、疑似距離を衛星時計のドリフトを考慮して修正する処理が行われます。 ところで、観測される擬似距離には本来なら受信機の時刻ドリフトも影響しますが、ここでは時刻同期の取れた受信機の疑似距離情報を使うので問題となりません。 推定された発信時刻を元に衛星の位置座標を計算し、さらに伝搬時間の分座標系を回転させて地球の自転を勘案します。 サンプルプログラムをそのまま実行すると、計算した衛星の位置座標とGPS受信機の計算した衛星座標と比較して殆ど一致することが確認されます。
実際の測位演算では、始めは疑似距離の計算なんてかなり誤差を含む訳ですから収束はより時間がかかるでしょう（たぶん）。 1回目の計算で求められた時刻の修正によって2回目はもう少し正確な衛星位置で計算することになります。 計算速度を短縮させるためには恐らく色々なアルゴリズムがあるのでしょう。
本節では、以上の様な測位に必要な衛星座標を求めるのではなくて、単純にあるエポックにおける衛星座標を求めてその精度を比較する事を目的とします。 そのため色々な補正処理を削除する必要があります。 取り敢えずは比較するための精密軌道暦を手に入れましょう。
サンプルとしてプログラムと同時に配布されているのは2002年4月25日の観測データとエフェメリス（一つのファイル内に$～で区別）です。 この日はGPS週番号1163の木曜日なので第4日です（日曜日が0）。という訳で、「http://igscb.jpl.nasa.gov/igscb/product/1163/」から「igs11634.sp3.Z」をダウンロードします。 このファイルはWindowsでは余り見かけない形式で圧縮されています。 この圧縮ファイルを解凍するソフトはLinux系ディストリビューションであれば標準装備ですがWindowsでこれを解凍するには何らかの解凍ソフトが必要です。 圧縮解凍ソフトでお勧めなのはLhaplusですが、インストール無しに使いたい場合であれば解凍レンジをお勧めします。 もし解凍ソフトをまだお持ちでなければvectorなどのサイトからダウンロードして下さい。 Lhaplusを使って解凍すると、“Z”の取れたフォルダが出来ています。 そのフォルダを開くと、解凍前と同じ名前のファイルが出来ています。 一見失敗したのかと思えますが、そのまま適当なエディタへドラッグ＆ドロップすることで閲覧することが出来ます。
ファイルのダウンロード画面は使用するブラウザによって異なります。 ファイル名の上でクリックするか、右クリックから保存を選べば保存が出来ます。
図2.4-1
図2.4-2
eph2ecef.cを実行すると時刻と衛星位置座標が表示されるようにソースコードを改造します。 出力させるのはGPS時刻・DOW（曜日）・時・分・PRN番号・座標X・Y・Zです。 デフォルトで出力されるメッセージを抑制するため、ヘッダファイルのDEBUG1を0に設定します。
（ ところで、衛星座標を計算する関数では擬似距離を与えるとそれだけ進む時間の分Z軸を回転させてサニャック効果を含めて座標を計算してしまいます。 便利な気もしますが、関数の汎用性を考えるとサニャック効果の分回転させるのは後処理でやった方が良い気がします。 ）
この改造したプログラムを実行すると、図2.4-6の様になりました。 図2.4-7が精密軌道暦です。 座標値の単位は精密軌道暦がkmで計算値の方はメートル（m）です。 ここで、14番衛星についてその位置座標を比較してみると何んと1.91m程度の誤差しかありません。 エフェメリスって、すごいんですね。 アンテナオフセットを入れたら、もっと高精度であることが分かるかもしれません。 ちなみに、受信機の出力していた衛星位置とは約270mずれていました。 これはサニャック効果と伝搬時間を考えたためです。 今、仮に伝搬時間が70ms，衛星軌道半径を26,000km，速度3.8km/s，衛星位置が赤道面とすれば、 伝搬時間の間に衛星は宇宙空間を266m進み、更に自転により見かけ上132m滑ります。 後は内積を取らないと分からないのですが、最大398m，最低134mものずれが生じることになります。
次の節ではさらにエフェメリスの精度について考察して行きます。
別ページにおいて、このプログラムを改造してエフェメリスの考察を続けます。
 図2.4-6 計算プログラムの実行結果 上段：サニャック効果も、エポック時刻からの伝搬時間分の引き算もやっていない結果 下段：デフォルトの計算結果 |
 図2.4-7 左と同時刻における精密軌道暦 |