2010/6/12更新
GPS単独測位プログラムの理解について
本サイトでは、GPS単独測位プログラムを理解するための
簡単な説明を行います。
3月30日の打ち合わせで話したように
ソースコードは、坂井さんのプログラムを利用します。
教科書は研究室に2、3冊ありますので参考にしてください。
ソースは電機大出版よりダウンロードしてください。
全体の流れは、以下の通りです
① 単独測位の概観
単独測位計算部分は、GPS受信機内部の処理でいうと、最後のほうになります。
単独測位に必要となる主なデータは、「疑似距離と衛星位置計算に必要なエフェメリス」です。
実際の受信機内部処理は、疑似距離やエフェメリスをきっちり計算するための信号処理が大部分で、
測位計算自体は、それらでてきた結果を利用して計算する部分になります。かといって
重要でないということはなく、ここをきっちりやらないと実際のアプリケーションで
利用される緯度、経度情報等の精度がおろそかになってしまいます。
② 単独測位に必要な生データ
単独測位に必要なデータは、上にも書きましたが、疑似距離とエフェメリスです。
③ その生データの入手方法
生データは、国土地理院のサイトより入手できます。ただし全て静止データとなります。
国土地理院→基準点・測地観測データ→電子基準点→提供サービスとすすんでください。
ここで、同意してデータが入手できるようになります。
てはじめに、2010年1月1日の千葉市川のデータを入手してください。
電子基準点検索より、日本、関東と進んで市川付近の赤い点をクリックし、電子基準点情報
にすすんでください。開始日と終了日を2010年1月1日として、観測データ閲覧をクリック。
ボタンにチェックをいれて、観測データ詳細一覧をクリックすると、ダウンロード画面になります。
2つのデータをダウンロードとして、解凍し保存しておいてください。
さらに、提供サービスのトップにもどり、左欄のデータの所在を知りたいときをクリックしてください。
右の日々の座標値をクリックしてください。日々の座標値「F3」の取得をクリックしてください。
さらに2010年のディレクトリへ。さきほど取得した市川基準局の番号は、93023です。その番号の
ファイルをクリックし保存してください。93023.10.posです。そのファイルを開いて、
2010年1月1日の部分をみて、地球中心のX,Y,Zの値と緯度経度高度の値を別の場所に
コピーしておいてください。この数値を、市川基準局アンテナの真値とします。
精度はよくわかりませんが、基本的に1cm程度で管理されていると思います。よって
単独測位の評価には十分利用できます。
4月12日までの課題を以下に書きます。
１）プログラムをダウンロードして、DOS画面で動作確認をしてください。
そのときに結果ファイルをきちんと出力してください。
pos2.exe ****.**n ****.**o > ****.csv
でよいと思います。
生データはサンプルについているデータと上記で取得したデータ両方必ずトライしてください。
２）次に自身の持つVisualStudioでプロジェクトを作り、そのプロジェクトにpos2.cを読み込んでください。
すぐにリビルドまでできるか確認してください。やり方がわからないときは必ず質問してください。
以下に簡単に示します。
新規作成で新しいプロジェクト、win32コンソールアプリケーション、場所と名前は自分でわかりやすいように。
アプリケーションの設定を選択し、コンソールアプリケーションと空のプロジェクトをチェックしてください。
これで完了です。次に今回作成したこのアプリケーションのファイルのある場所に
pos2.cをコピーします。そして、VisualStudioのプロジェクトから既存項目の追加を選び、その
pos2.cを選びます。これでプロジェクトに追加されます。
３）次に、このままでは実行できないので、main文の中を少し
書き換えてみてください（main文は一番下のようにあります）。main文は
必ず全て理解するようにしてください。理解できれば修正も簡単になります。
基本的に2つのファイルを読み込むように
設定すればOKです。現状は、コマンドラインでオプションと2つのファイルを
読み込む設定になっています。これを変更して2つのファイルをmain文に書いて
読み込むような形にしてみてください。
そうすると、DOS画面と同じように実行できます。
この実行までを12日までに確認してください。
なお、実行の際に2つのファイルはpos2.cと同じ場所においたほうがよいでしょう。
ディレクトリを指定する必要がありません。
4月16日までの課題を以下に書きます。
１）まず上記の３）までの課題を確実に行ってください。main文を少し
変更することで可能になります。
２）次に、これまでコマンドラインに出力されていた測位結果を
ファイルに出力するように変更しましょう。
具体的には、compute_positionの下のほうで測位結果が出力されているので
そこを書き換えます。一番下のほうにelse文があり、受信機位置等が
出力されている箇所があります。この部分をコメントアウトして
新たにfprintfで全く同じように出力させます。
fprintfの使用方法は、C言語の教科書等を参考にしてください。
簡単にできるはずです。
３）ファイルに書き出せるかどうか、存在する三鷹と市川のデータで確認してください。
次に、市川のデータよりアンテナの真値がわかっているはずなので、
それをもとに、地球中心でX,Y,Z方向の標準偏差を計算してください。
それを16日に提出してください。
4月27日までの課題を以下に書きます。
測位計算を行う以下の関数で、仰角による使用衛星判別を行っています。まずその部分を探してください。
次に、実際に仰角で判定を行う際の閾値がありますので、その変数を見つけてください。
さらにその閾値を30度とした場合と0度とした場合の
測位結果を比較してください。single_pos.csvに出力される
可視衛星数の1日変化を必ずグラフにして観察するようにしてください。
なお、仰角による使用衛星判別を行っている関数名は以下です。
static int _compute_position(wtime wt,double *psr1,bool detail,
double sol[MAX_M],double cov[MAX_M][MAX_M],double *dpsr,
double *dpsr1,double *el,double *az)
④ 生データを読み込む部分
ここは観測データと航法メッセージの読み込み部分なので、
どのように読み込んでいるかだけチェックしておいてください。
⑤ 衛星位置計算方法
5月17日までの課題を以下に書きます。
衛星位置計算を行う関数と衛星のクロックを補正する関数内で下記のことをためしてください。
まず、衛星クロックの補正する関数「satellite_clock」の関数内に「衛星時計の補正量を計算」
とあります。この部分は、衛星クロックがいくら原子時計を積んでいても
1ns以内の正確さがないことを示しています。
そこで、この関数のdtを0にした場合とそのまま計算した場合で
どの程度測位精度に影響があるかをチェックしてください。
次に、衛星位置計算する関数「satellite_position」の関数内に「補正係数を適用する」
とあります。この部分は、アルマナックとエフェメリスの大きな差を生んでいます。
おおざっぱですが、このd_uk,d_rk,d_ikの3つの変数を0にすると補正されないことになり
アルマナック相当の精度になってしまいます。
この部分もまた、どの程度測位精度に影響があるかをチェックしてください。
⑥ 衛星の仰角・方位角計算方法
実際に衛星の仰角や方位角を計算するには、衛星位置を算出するときに利用している
地球中心座標系では無理。地球中心座標系は、地球中心が（0,0,0）で、
基準となるユーザ位置と衛星位置もこの原点に基づいています。
単純に考えれば、その座標のまま仰角や方位角を求めることは困難です。
では、どのようにすればユーザである自分のアンテナに対する衛星の仰角と方位角を
求めることができるかというと、自分のアンテナを中心にした座標系を
作ればよいことがわかる（アンテナを頂点とする地球表面で切る）。自分のアンテナを原点（0，0，0）とし、
北方向をy軸、東方向をx軸、天頂方向をz軸とする座標系を想定すれば、
自分のアンテナに対する衛星の仰角と方位角は簡単な計算より算出することが
可能となる。そのような座標系を求めるために、地球中心だった座標系を
変換する必要がある。その手順は以下の通りです。例えば、日本を東経35度、北緯35度、地面をアンテナ位置とします。
（地球中心座標系では、経度0方向をx軸、地軸方向をz軸、経度90度方向をy軸としています）
１）地球中心座標系をまず自分のアンテナ位置まで平行移動します。
２）次に、地軸方向を向いている軸が、天頂方向を向くように（90－緯度）分だけ
回転させます。緯度が35度の場合55度赤道側に倒すイメージ。
３）次に、経度0度方向を向いているx軸を地面上で東方向に向くよう
右に経度分（135度）だけ回転させます。こうすることにより、
自動的に経度90度方向を向いていたy軸が北方向を向きます。
上記の3つの移動及び回転を地球中心座標系の衛星位置（x,y,z）に対して
行うと、自動的に自分のアンテナ位置を中心とした座標系に
それらの（x,y,z）が変換されます。そして、新しい（x,y,z）を利用して
仰角や方位角を求めることになります。
実際に坂井さんのプログラムでは、posenu_xyz_to_enu(posxyz pos, posxyz base)
という関数で求められています。この関数内での算出方法は上記とは
やや異なりますが、考え方はいっしょです。
pos.x,pos.y,pos.zは構造体で、衛星の地球中心座標位置が入っています。
base.x,base.y,base.zは構造体で、ユーザの地球中心座標位置が入っています。
この関数では、オリジナルの相対的なベクトル（上記のx,y,z軸方向の差分）を変換しています。
回転行列に関しては、みなさんの教科書のP76,77に記述があります。
左回りなのか右回りなのかに注意してください。
坂井さんのプログラムによる変換方法を自分で確認しておいてください。
5月24日までの課題は以下の通りです。
１）上記のプログラムの中身を理解した後、実際に下記の座標について
仰角と方位角を手の計算で求めてください。提出は24日にレポートでOKです。
x,y,z軸の3次元の座標系があるとし、原点を（0,0,0）とします。
この座標系でユーザは（0,1,1）に位置し、衛星は（-2,2,2）に位置するとします。
このとき、ユーザに対する衛星の仰角と方位角を求めてください。
ユーザ位置の緯度経度高度は、北緯45度、東経90度、高度0mとします。
レポートは手書きが望ましいですが、電子ファイルでもOKです。
x,y,z軸の3次元の座標系のイメージ図を以下に貼っておきます。
元の座標系イメージ図
⑦ 電離層や対流圏遅延量の推定
電離層遅延量はiono_correctionで計算されています。
対流圏遅延量はtropo_correctionで計算されています。
それぞれ、教科書に書いてあったアルゴリズムと同様に
計算されていることを確認してください。教科書でいうと、
電離層は、P105から106付近、対流圏は、P189付近になります。
ただし、実際に使用されている計算式を別の教科書でも見ておいてください。
研究室にあります（精説GPS）。
5月31日までの課題は以下の通り。
１）電離層遅延量を0にした場合の測位結果を見てください。
電離層遅延推定ありとなしの場合で、水平及び高度方向の測位精度を比較してください。
２）対流圏遅延量を0にした場合の測位結果を見てください。
対流圏遅延推定ありとなしの場合で、水平及び高度方向の測位精度を比較してください。
３）電離層及び対流圏遅延量双方を0にした場合で、水平及び高度方向の測位精度を比較してください。
４）時系列の電離層遅延量と対流圏遅延量を時刻、仰角とともに
ファイルに書き出してください。そして、
仰角に応じてどのように遅延量が変化するかグラフに示し確認してください。
31日に、上記の精度比較をまとめたレポート1枚を提出してください。
測位精度の比較は、真値からのずれ平均と標準偏差でOKです。
（地球中心座標系のx,y,z方向の結果でOKです。可能であれば緯度、経度、高度方向の測地座標系に直してください）
実際にファイルに書き出す部分は慣れないかもしれませんが、
重要なので書き出し方法がわからない場合は、上記を復習するか
久保まで質問してください。
24日の課題の注意事項ですが。教科書P76の回転行列は軸方向を見て
時計周りが正となっていると思います。時計周りを正とするか負とするかに注意してください。
(回転行列については、軸周りの符号によってsinの前の符号が変わるのかと思います)
⑧ 測位に使用する衛星を決定する部分
こちらは、最低仰角や最低信号強度を決定して、測位計算に使用する
衛星を決定する部分です。他にも様々なファクターがあります。
信号処理のTracking部をパスした衛星については
全て擬似距離、搬送波位相等（最低でも擬似距離）が出力されています。
しかしながら、全ての衛星を利用すればそれでいいのかということではなく、
例えば、仰角が0度付近の衛星はやはりマルチパス誤差が
それなりに大きいです。例えば、最低仰角15度とし、8機の使用衛星で水平で1－2mの測位精度のときに
仰角0度付近の衛星2機を測位計算に入れて、10機とすると数m逸脱することもあると思います。
必ずしも衛星配置が全てではないということです。もちろん
衛星配置が良いにこしたことはありませんが、あくまでも測定誤差とのバランスが重要です。
よって、アプリケーションによってどのように閾値を決定するかが
重要になります。例えば、飛行機のナビゲーションに利用する場合は、
この部分が重要になることはわかりますね。
もちろん受信機単独のチェックでは信頼性等は担保されないため、
他のシステムによる2重、3重のチェック機能がついてくると思います。
通常の測位計算では、マスク角10－15度、最低信号強度30dB-Hz程度でOKです。
⑨ 実際の単独測位
実際の単独測位は、compute_positionで行われています。compute_positionの関数を
上からみていくと、まずエフェメリスをセットして、解の初期化を行っています。
次に3回荒い計算を行っています。実際の計算関数は「_compute_position」です。
赤いbreakポイントで_compute_position関数の上のほう（暫定のユーザ位置）で停止させると、
地球中心の解（x,y,z）が（0,0,0）になっていると思います。確認してください。
これは、最小二乗法の初期値を地球中心の原点としていることになります。
次に_compute_position関数の一番下で停止（break）させてください。
1回目の計算で解（x,y,z,dt）が（-4533515.53,3924881.85,4325320.95,1117203.27）
となっていることを確認してください。単位はmです。これらの値はまだ日本付近にも
きていません。最小二乗法は、与えられた複数の衛星位置と、複数の距離情報から
最もらしい位置を求めるために利用しています。一番落ち着く場所を探しているとも
いえます。逐次法のため、収束するまでに何回か計算する必要があります。
ここで、上記で述べた、与えられた衛星位置と与えられた距離情報がもし完璧な
真値だった場合は、最小二乗法で計算しても必ず真値の1点に求まります。
2回目の計算で解（x,y,z,dt）が（-3964303.56,3381864.26,3715326.12,30897.73）
となっていることを確認してください。なんか実際のアンテナ位置に近づいてきたように思います。
初期値を地球中心の原点というとんでもない位置にしたとしても、2回の計算で
ここまで近づいてきます。
3回目の計算で解（x,y,z,dt）が（-3947783.64,3364421.28,3699449.57,27.07）
となります。
4回目以降は、_compute_position関数の精密計算という部分に入ることを確認してください。
これは数cmの精密という意味ではなく、ここから電離層や対流圏も推定して
距離情報をより正確したもので計算することになります。確認してください。
4回目の計算で解（x,y,z,dt）が（-3947762.76,3364400.16,3699429.37,-14.15）
となります。
5回目の計算で解（x,y,z,dt）が（-3947762.74,3364400.14,3699429.35,-14.19）
となります。もう4回目と比較しても数cmしか動いていません。
だいぶ収束してきました。
6回目の計算で解（x,y,z,dt）が（-3947762.74,3364400.14,3699429.35,-14.19）
となります。4回目と比較してcmの単位ではもう動いていません。
7，8回目も同じようになることを確認してください。最小二乗法で求められる解の
前回との誤差がある一定値以下になると終わるようにしていると思います。
（私は確認していません）
以上ですが、この部分は単独測位を知る上で重要な部分ですので、
しっかり理解しておいてください。
21日までの課題は以下の通りです。
今回の計算で初期値の地球中心から回を重ねるごとに、推定位置が
どのような場所に推移していったかをきっちり確認してください。
できればある国のある地域の地下方向何m（または標高何m）という形で。
地球中心のx,y,zを緯度、経度、高度に直せばわかるかと思います。
次に初期値を（0,0,0）にせずに、どこか異なる場所にしたときの
解の推移を上記と同じようにレポートにして提出してください。
回数は5回程度まででOKです。
不明な点がある場合は必ず連絡ください。
⑩ 受信機時計誤差の取り扱い
⑪ 様々なデータを解析してみよう
by 久保