電子基準点 小松(0255) その(2)

実家から帰宅。運転疲れた。

「小松」にせっかく行ってきたので、マルチパスを調べてみることにした。
MPを計算するプログラムも書いてみたことだし、チェックを兼ねて。

ということで、実家近くの「小松」（上）と「白峰」（下）を比べてみたところ。
局座標はアプリオリでズルしている。

１．「小松」vs「白峰」　（スカイプロット：小松（上）、白峰（下））
国土地理院のページにある「スカイプロット表示」を 真似たもの。

２．「小松」vs「白峰」　（仰角横軸：小松（上）、白峰（中）、重ねたもの（下））

スカイプロットを見ると「小松」が圧倒的に観測条件は良いのだけど、仰角で見てみると、高仰角では「白峰」の方が良いような。

「小松」が海沿い（海抜0m近い）にあるのに対して、「白峰」は500mくらいのところにあるため、大気遅延誤差が軽減されているんだろうか。
（「小松」のMP1で仰角7度あたりにトビがあるのは、平均を取る際に適当に仰角10度以下の値を切ったせい？？）

とにかく、「小松」は良い観測点のようで良かった。地元民として。

さて、 スカイプロット表示は以外にややこしかったのでgnuplotスクリプトを貼っておく。

set polar
set size square
unset border
set format x ""; set format y ""

set angles degrees
set datafile commentschars "#"
set datafile separator ","
set xrange[-90:90]
set yrange[-90:90]
set cbrange[-6:6]
set grid polar 30
set xtics axis 15
set ytics axis 15

# labels
set label "15" at 15*cos(120),15*sin(120)
set label "30" at 30*cos(120),30*sin(120)
set label "45" at 45*cos(120),45*sin(120)
set label "60" at 60*cos(120),60*sin(120)
set label "75" at 75*cos(120),75*sin(120)
set label "90" at 90*cos(120),90*sin(120)

set palette rgbformulae 22,13,-31

set multiplot title "Code multipath" layout 1,2
plot "tmp.dat" using (90+$1):(90-$2):3 pt 7 ps 0.3 lc palette title "MP1"
plot "tmp.dat" using (90+$1):(90-$2):4 pt 7 ps 0.3 lc palette title "MP2"
unset multiplot

tmp.datの中身はこんな感じで方位角、仰角、MP1、MP2が並んだもの。

 75.277634,  11.737041,  0.329021,  -1.379482
 75.450305,  11.595049,  0.588193,  0.520832
 75.622963,  11.453226,  -0.309160,  0.994863
 75.795605,  11.311571,  0.180602,  -0.183782

電子基準点 小松(0255)

お盆で帰省中なので、家の近くの電子基準点の周辺状況を確認した。
なかなか良いところに設置されていましたが、基礎とのつなぎ目がすこし割れていた。

後ろに松林（東側） がありますが、受信障害になるレベルではなさそうですね。
いいところですな、石川県小松市。

ちなみに場所はRTKLIBで計算してgoogle検索しました。（またRTKLIBかと）

衛星のephemerisの扱い

単独測位をやってみようか、と思ってふと気づいたのだけども、受信機が記録している時刻は、衛星の電波発信時刻と異なる（いまさら気づく）。

いくら光速で信号が届いたからと言って、受信時刻と同じ時刻なわけはない。

さあ、普通はどのように扱うのかな、と思って先日買った教科書「GNSS」を開く。

えー、
"受信機まで0.06秒くらいかかります。衛星は1km/sで動いているので60mに相当します。"

おっと、だめじゃないか。QZSだったりすると軌道がでかいからさらに大きくなるぞ。しかも教科書に補正方法書いてないじゃん。

ということで、ephemerisを座標値に変換する場合は、多分次のようにする必要がある。

暦計算の時刻＝受信時刻 −（擬似距離÷光速）＝発信時刻

でもこの程度でいいのかな？直接軌道誤差に効いてきてしまう。

あと、暦の中の時計誤差は、こうやって出した時刻に対しておこなって、その時刻で算出すればいいか。つまり二回時計誤差の計算をすることになるかな。

IS-GPSにきちんと書いてあるかもしれないけど読む気力が無い・・・

まあとにかく受信時刻で暦を使っちゃだめだよ、と今更気づいた、という話。

MP線形結合の計算（２）

昨日MP線形結合の計算をしてみたけども、出力してみるとやけに大きな値が出る衛星がある。

なんだろうな、と思ったらGLONASSだった。GLONASSのFDMA波長にきちんと衛星毎に合わせる必要があるのだった。

少し計算を分けてみたのでメモしておく。
例によってrtklibに便利なsatwavelenという関数があることに気づいたので、早速使ってみる。

/* compute mp linear combination for GPS and QZS */
static int mp_gps(const double P1, const double P2, const double L1, const double L2, double *mp1, double *mp2) {
    double a=FREQ1*FREQ1/(FREQ2*FREQ2);
    double lam1=CLIGHT/FREQ1;
    double lam2=CLIGHT/FREQ2;

    *mp1 = P1 - (1+2/(a-1))*L1*lam1 + (2/(a-1))*L2*lam2;
    *mp2 = P2 - (2*a/(a-1))*L1*lam1 + (2*a/(a-1)-1)*L2*lam2;

    return 0;
}

/* compute mp linear combination for GLONASS FDMA */
static int mp_glo(int satnum, const nav_t *nav, const double P1, const double P2, const double L1, const double L2, double *mp1, double *mp2) {
    double a, lam1, lam2;

    lam1 = satwavelen(satnum, 0, nav);
    lam2 = satwavelen(satnum, 1, nav);
    a = lam2*lam2/(lam1*lam1); /* f1^2/f2^2 */

    *mp1 = P1 - (1+2/(a-1))*L1*lam1 + (2/(a-1))*L2*lam2;
    *mp2 = P2 - (2*a/(a-1))*L1*lam1 + (2*a/(a-1)-1)*L2*lam2;

    return 0; 
}

aは、f1^2/f2^2だけど、波長だと分母分子が入れ替わります。
satwavelenがすごく便利でした。このコード例では毎回呼び出すのでちょっと無駄な処理が多いんですが。

出力：
----- epoch start -----
sat num: 17, timetag, PRN, MP1, MP2
 0: 2012/07/09 00:00:00   9       71.8339      100.7876
 1: 2012/07/09 00:00:00  41     -242.7161     -297.4944
 2: 2012/07/09 00:00:00  14      137.7746      184.8548
 3: 2012/07/09 00:00:00  40     -242.3190     -297.2934
 4: 2012/07/09 00:00:00  15      119.1687      160.6320
 5: 2012/07/09 00:00:00  42      377.8123      480.3067
 6: 2012/07/09 00:00:00  31       61.5032       91.2867
 7: 2012/07/09 00:00:00  12       35.9511       58.4904
 8: 2012/07/09 00:00:00  33      183.4078      236.6984
 9: 2012/07/09 00:00:00  22       28.5107       45.0275
10: 2012/07/09 00:00:00  25       21.7292       46.4806
11: 2012/07/09 00:00:00  18       77.3021      105.7750
12: 2012/07/09 00:00:00  43     -527.0718     -645.2571
13: 2012/07/09 00:00:00  51       34.2777       53.0246
14: 2012/07/09 00:00:00  52     -118.9257     -130.6246
15: 2012/07/09 00:00:00  27       98.0804      133.2339
16: 2012/07/09 00:00:00  39       11.8226       19.8136

とりあえずまともそうな値が出たのでよしとする。多分合ってるだろう。
あとは単独測位をして仰角を横軸にしてノイズレベルを調べてみようかな。

MP線形結合の計算

MP線形結合、と呼ぶのかどうかわからないけど、teqcのマニュアルを見ると載っているクオリティチェック用の線形結合。

要は、電離層遅延と対流圏遅延の両方を差っ引いて、マルチパスノイズの変化量だけにしたもの。

計算するだけなら簡単。ということでメモしておく。
例によって、rtklibを使ったもの。

/* print code multipath for L1, L2 */
int printmp(const obs_t *obs) {
    int i=0, j, nu;
    char tmstr[20];
    double mp1, mp2;
    double a=FREQ1*FREQ1/(FREQ2*FREQ2);  /* alpha = (f1/f2)^2 */

    while(in) {
        /* seek index of the next epoch */
        nu = nextobsf(obs, &i, 0);
        printf("----- epoch start -----\n");
        printf("sat num: %d, timetag, PRN, MP1, MP2\n", nu);

        for(j=i;jdata[j].time, tmstr, 0);
            mp1 = obs->data[j].P[0] - (1+2/(a-1))*obs->data[j].L[0] + (2/(a-1))*obs->data[j].L[1];
            mp2 = obs->data[j].P[1] - (2*a/(a-1))*obs->data[j].L[0] + (2*a/(a-1)-1)*obs->data[j].L[1];

            printf("%2d: %s %3d %13.4f %13.4f\n", j-i, tmstr, obs->data[j].sat, mp1, mp2);
        }
        i += nu;
    }


    return 0;
}

よく見ると、メートルに直し忘れているけれどそのままにしておく。

ubuntu12.04でのgpscalの不具合(でもないけど)

国土地理院のftpサーバで、"gpscal"という便利なスクリプトが公開されている。

gpscal:
ftp://terras.gsi.go.jp/software/gpscal/

これは、linuxのcalコマンドの出力をgps week等を表示するよう整形するスクリプト。
とても便利なんですが、どうやらスクリプト内で使用しているcalコマンドの仕様がubuntu12.04では変更になってしまったらしく、現在日の表示がうまくいかない。

例えば、週の頭に実行してしまったりすると、丸々週が抜ける。

              August 2012         
Week   Sun Mon Tue Wed Thu Fri Sat
1699                 1   2   3   4
                   214 215 216 217
1700    12  13  14  15  16  17  18
       225 226 227 228 229 230 231
1701    19  20  21  22  23  24  25
       232 233 234 235 236 237 238
1702    26  27  28  29  30  31    
       239 240 241 242 243 244

これは、calコマンドが現在日をハイライトしているためにうまくawkで処理できなくなっているためなので、次のように"gpscal"内の"cal"呼び出しオプションに"-h"を追加すればOK.

cal -h $mm $yyyy |awk 'BEGIN {doy='"$doy"'; week='"$week"';\

しかしcalコマンドも派生版がたくさんあるようなので、環境毎に対応するしかないのかな。今の美しいスクリプトを汚くするのも嫌だし。