対流圏遅延補正

今日は対流圏遅延補正のHopfieldモデルと、Saastamoinenモデルについて少し考えてみた。
どっちも良く聞くのだけども、どのように差が出るのかとか、あまり考えたことがない。

ということで、前買った教科書「GNSS」に書いてあるモデルの式を試しに計算してみた。何が疲れたって、Saastamoinenの式の高度補正のテーブルを打ち込むのが疲れた。

さて、海抜0kmから5kmまでを計算して比べてみる。

（実はテーブルの補正値が5kmまでしか載っていないからなんだけども）

計算条件：
大気圧　P=1013.25 (hPa)
分圧　P(dry)=1000.00 (hPa), P(wet)=10.00 (hPa)
気温　T=300 (K)
分圧の和がおかしいとかはおいておき。

(Hopfiled model)

 (Saastamoinen model)

比べてみると、0kmの標高ではほとんど値は同じ。
Hopfieldモデルでは標高が上がるにつれて遅延量が小さくなる。
これは、通過する大気の厚さが減るからで、イメージ通り。

だけど、Saastamoinenのモデルで値がほとんど変わらないのはなぜ？？
少し拡大してみると、

 

一応、値はちゃんと変わっているんだけど、Saastamoinenの式への補正量があまり効いていない。実際見積ると5km変わっても大して影響しない。

仮に計算が合っているとすれば、大気圧がこのモデルの支配的なパラメータになっているからじゃないだろうか、と想像。標高というより、大気圧を水蒸気量に換算してるんじゃないかと思う。元論文を読んでみたい・・・。

結局、まともな値にするためには大気圧を標高によって計算しなおしてあげないといけなさそうな気がする。

さて、標高0kmだとほとんど差がないんだけど、実際はどちらを使うのがいいのだろう。


あと、参考に、実測値を探していたら海外の文献がgoogleで引っかかった。
 http://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0CCMQFjAA&url=http%3A%2F%2Fwww.fig.net%2Fpub%2Ffig2006%2Fpapers%2Fps05_04%2Fps05_04_09_katsougiannopoulos_0725.pdf&ei=hsCDUIPGHevsmAXlmoGQBw&usg=AFQjCNHBxdIrHcEwMCzxPGp3MvAPClMp2w&sig2=DKD-zTzAFhdO1qIKa8f9HQ

比べると大体同じなので、プログラムに大きなバグがあるわけではなさそう。ほっとした。

prettyprintできれいなコードを

この前、bloggerにきれいなコードを貼るjavascriptを導入した。
"prettyprint"というやつでgoogleから無料でダウンロードして自前のサーバに置いてもいいし、googleのサーバに直接見に行くこともできる。
bloggerレイアウトのヘッダ部分に指定すればOK。
http://google-code-prettify.googlecode.com/svn/trunk/src/prettify.css
http://google-code-prettify.googlecode.com/svn/trunk/src/prettify.js

ただ、いつも書き始めると忘れているので使い方をメモしておく。

<pre class="prettyprint" style="font-size: medium;"><code>

こんな風にclassを指定するだけ。
すばらしい。

電離層遅延（放送暦）

今日はGPS放送暦の電離層遅延について少し調べた。
IS-GPS-200Eを読んでみただけだけど、混乱したのは、角度の単位が式で異なること。間違っているかもしれないけれども、また混乱しそうなので、一応単位をメモしておく。

F=1+16*(0.53-EL)^3 (EL:semi-circle)
phi=0.0137/(EL+0.11)-0.022 (EL:semi-circle)
cos(AZ), sin(AZ) (AZ:radian)

もちろん三角関数の中身がsemi-circleになるはずはないのだけど、Fなんかはどっちだ?となってしまった。

さて、IS-GPS-200Eにのっているモデルはklobucharモデルというやつなのだけど、元論文を読んだことはない。だけど、式を見ると、「何時にこの座標でこれだけの遅延量がある」という記述というより、「丸っこい遅延の塊が移動する様子」が記述されているように見える。なるほど・・・、電離層遅延は伝搬していくということなのかな?

勉強はこれでやったことにして、国土地理院からダウンロードしたRINEXファイルを確認してみると・・・、なんと、電離層遅延パラメータが載っていないではないですか！RINEXファイルフォーマットを読むと、オプション項目なんですね。

ふーむ。IGS点のデータを持ってくるか・・・。GEONETは定点観測なので、相対測位が前提ということなのかなあ。

まあ、とりあえずIS-GPS-200Eの該当部分を読んでプログラムにしてみたので、まだ動くかはわからないけど、超ベータ版ということで自分のために載せておくことにした。

double brdion(const double az, const double el, const double phi_u, const double rambda_u, const gtime_t time, const double *ion_model) {
    double dt;  /* ionospheric delay (sec) **not correction value** */

    double phi_i0, phi_i, rambda_i, phi_m, psi;
    double amp, per, f, x;
    gtime_t tl;
    const double R2SC=R2D/180.0;
    const double D2SC=1/180.0;
    const double SC2R=PI;
    int i;

    amp=per=0.0;

    /* estimate phi */
    psi = 0.0137/(el*R2SC) - 0.022;
    phi_i0 = phi_u*D2SC + psi*cos(az); /* az(rad) */
    
    if (fabs(phi_i0) <= 0.416) {
        phi_i = phi_i0;
    }
    else if (phi_i0 > 0.416) {
        phi_i = 0.416;
    }
    else {
        phi_i = -0.416;
    }

    rambda_i = rambda_u*D2SC + (psi*sin(az)/cos(phi_i*SC2R));
    phi_m = phi_i + 0.064*cos(rambda_i*SC2R - 1.617);

    f = 1.0 + 16.0*(0.53-el*R2SC)*(0.53-el*R2SC)*(0.53-el*R2SC);

    /* ion_model: iono model parameters {a0,a1,a2,a3,b0,b1,b2,b3} */
    for (i=0;i<4;i++) {
        amp += ion_model[i]*pow(phi_m,i);
        per += ion_model[4+i]*pow(phi_m,i);
    }
    if (amp < 0.0) {
        amp = 0.0;
    }
    if (per < 72000.0) {
        per = 72000.0;
    }
    tl = timeadd(time, 12.0*3600.0*rambda_i);
    x = 2.0*PI*(tl.time+tl.sec-14.0*3600.0)/per;

    if (fabs(x) <= 1.57) {
        dt = f*5.0e-9;
    }
    else {
        dt = f*(5.0e-9 + amp*(1-x*x/2.0+x*x*x*x/24.0));
    }

    return dt;
}

単独測位を試す

今日、昼にラーメンを食いに行こうと思ったら、前から行こうと思っていた店が休み、２番目に行った店は開店３０分前にもかかわらず行列が・・・。ラーメンは諦め、結局行きつけの店でカツ丼を食べました。女将さんからは「いつものね？」と言われた。相変わらず、おいしかったです。満足。


で、帰ると微妙な時間だったので、今日は単独測位プログラムを書いてみた。
観測方程式はこんな感じ。

 f(x,y,z,dt)=sqrt(dx^2+dy^2+dz^2)+C*dt

fは擬似距離から時計誤差とかを処理したもの。
x,y,zが線形じゃないので、これをテイラー展開して、非線形最小二乗法にする。
数式はノートに書いてみたら２ページになった上に、計算ミスもあったので今度メモっておきたいなあ。

さて、使ったデータは、GEONET「3018」観測点で、答えは(-3.9695074060E+06,3.3278898661E+06,3.7090383081E+06)　(2012/8/16のF3解)。

で、初期値を(0,0,0,0)から計算してみると、

satnum, x(sat), y(sat), z(sat), dt(sat)
SATPOS:  32,    8075372.3045,   22569776.5255,   12147617.4769, -0.0004770960
SATPOS:  84,  -23318319.2847,   18586114.8063,  -25160339.5912,  0.0000130641
SATPOS:  30,   -3332934.0086,   24125709.4567,   10248378.3800, -0.0003655502
SATPOS:  31,   -9053121.7708,   11965903.1817,   22118916.5717,  0.0002678009
SATPOS:  29,  -21467145.5567,     171477.0652,   15665776.8263,  0.0003446830
SATPOS:  16,     563918.4785,   26380546.3717,    -191508.9968, -0.0002438091
SATPOS:  25,  -15055889.6059,   -4018044.5479,   21474952.5222,  0.0000398373
SATPOS:  14,  -21625311.0310,   14935548.3151,    3395049.7157,  0.0002101171
itrnum, time, x(sta), y(sta), z(sta), dt(sta)
POS(000): 2012/08/16 00:00:00,-4755659.5459, 3686366.5401, 4073518.8387,0.003110526
POS(001): 2012/08/16 00:00:00,-4152400.3866, 3342603.6460, 3754318.7889,0.003189483
POS(002): 2012/08/16 00:00:00,-4143561.5152, 3338589.6516, 3752314.5547,0.003225822
POS(003): 2012/08/16 00:00:00,-4143559.3173, 3338587.5856, 3752310.6861,0.003262148
POS(004): 2012/08/16 00:00:00,-4143559.3125, 3338587.5793, 3752310.6798,0.003298474
POS(005): 2012/08/16 00:00:00,-4143559.3125, 3338587.5793, 3752310.6798,0.003334800

となった。５回くらいiterationかければ収束するのかぁ。

しかし座標値が２００kmとかずれている・・・。一体何なんだ。それなりに正解に近いというのがいやらしい。

ということで試行錯誤の末。
・・・衛星時計誤差を補正していないためだった。忘れていたよ・・・。
しかし、たった1ミリ秒以下じゃん、とたかをくくっていたけども、こんなに効くんだなぁ。へー。

ということで、こちらが修正後。

（衛星時計誤差を考慮）
satnum, x(sat), y(sat), z(sat), dt(sat)
SATPOS:  32,    8075372.3045,   22569776.5255,   12147617.4769, -0.0004770960
SATPOS:  84,  -23318319.2847,   18586114.8063,  -25160339.5912,  0.0000130641
SATPOS:  30,   -3332934.0086,   24125709.4567,   10248378.3800, -0.0003655502
SATPOS:  31,   -9053121.7708,   11965903.1817,   22118916.5717,  0.0002678009
SATPOS:  29,  -21467145.5567,     171477.0652,   15665776.8263,  0.0003446830
SATPOS:  16,     563918.4785,   26380546.3717,    -191508.9968, -0.0002438091
SATPOS:  25,  -15055889.6059,   -4018044.5479,   21474952.5222,  0.0000398373
SATPOS:  14,  -21625311.0310,   14935548.3151,    3395049.7157,  0.0002101171
itrnum, time, x(sta), y(sta), z(sta), dt(sta)
POS(000): 2012/08/16 00:00:00,-4546577.0224, 3674293.9179, 4025366.1831,0.002756903
POS(001): 2012/08/16 00:00:00,-3977504.4200, 3331504.2787, 3710605.1526,-0.000159663
POS(002): 2012/08/16 00:00:00,-3969535.1170, 3327893.6137, 3709067.7009,-0.000198327
POS(003): 2012/08/16 00:00:00,-3969533.8351, 3327893.1039, 3709067.8222,-0.000198333
POS(004): 2012/08/16 00:00:00,-3969533.8351, 3327893.1039, 3709067.8222,-0.000198333

ということで、大体数十メートルくらいの精度になったみたいだ。
まだ相対論効果、電離層遅延、対流圏遅延を考慮していないけど、これくらいにはなるんだね。
次はこれらを入れたらどう変わるかを見てみようかな。

今日は、単独測位プログラム程度とは言ってもそれなりに計算しなきゃいけないんだな、ということが分かった。相対測位とか面倒くさそうだな、こりゃ。