ついにRTKLIB2.4.2ベータ公開！！

待ちに待ったRTKLIB2.4.2（ベータ）がいつの間にか公開されているのに気づいた！！！

http://www.rtklib.com/
（なぜか更新日「2013 1Q」？？だけど）

まだリリースノートは出ていないのだけど、ざっとソースを確認してみる。
どうやら、相対測位のambiguity resolutionに「WL-NL」と「TCAR」が加わったようだ。
あと、PPP-ARが追加。こ、これは・・・。

PPP-ARのアルゴリズムは全く勉強してないので、さっそくソースから勉強させてもらおう。

「WL-NL」は一体どんな結果になるんだろうか。長距離までいけるんだろうか。
「TCAR」はまだL5があまり使えないので実際に試すことはできるのかな？

あれ、そういえばfix-and-holdと一緒に指定できないけど、この２つは平均値を丸めてambiguityを求めるからfix-and-holdで拘束をかけれないということ？うーん、うまいこと共存できないんだろうか。

ああ、しかしこんな時間（03:20）に気づくなんて！ のんびり酒飲んでニコ動見てる場合じゃなかったよ。だけど、呼び出し方だけとりあえずメモっておく。

ソースを読んだ感じでは、コンフィグファイルに次のように書きこめばambiguity resolutionの方法が切り替わるはず。

pos2-armode = 4 # widelane-narrowlane
pos2-armode = 5 # TCAR

明日（すでに今日）は起きたら早速試してみたい。
うおお、楽しみだー。

（追記）
電離層読み込み機能も追加されているみたい。
PPP-ARに使うのかな。

対流圏遅延補正

今日は対流圏遅延補正の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

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

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

逆行列の計算(fortran)

今日は逆行列の計算をしてみたのでメモ。

ソルバはLAPACKのDGETRF,DGETRIで、fortran90。

! --------------------------------
! calculate inverse of matrix A(m,n)
! --------------------------------
function inv_lapack(A) result (AI)
    implicit none

    real(8), intent(in) :: A(:,:)
    real(8) :: AI(size(A,1),size(A,2))
    integer, parameter :: NB = 64
    integer :: info
    integer :: m, n

    ! for lapack
    integer :: LWORK
    integer, allocatable :: IPIV(:) ! dimension N
    real(8), allocatable :: WORK(:) ! dimension LWORK

    ! get size of A(m,n)
    m = size(A,1)
    n = size(A,2)

    LWORK = n*NB
    allocate(IPIV(n))
    allocate(WORK(LWORK))

    AI = A
    call DGETRF(m, n, AI, m, IPIV(1:min(m,n)), info)  ! factorize
    call DGETRI(n, AI, m, IPIV, WORK, LWORK, info)  ! inverse

    deallocate(IPIV)
    deallocate(WORK)

    if (INFO==0) then
        !write(*,*) "succeeded."
    endif

    return
end function

! --------------------------------
! show matrix A(m,n)
! --------------------------------
subroutine show_dmat(A,fmtf)
    implicit none
    real(8), intent(in) :: A(:,:)
    character(len=*), intent(in), optional :: fmtf
    character(len=25) :: fmtc
    integer :: m, n, i, j
    m = size(A,1); n = size(A,2)

    if (present(fmtf)) then
        fmtc = fmtf
    else
        fmtc = '(f13.4," ")'
    endif
    do i=1,m
        do j=1,n
            write(*,fmtc, advance='no') A(i,j)
        enddo
        write(*,*) ""
    enddo
    return
end subroutine show_dmat

呼び出しはこんな感じ。ついでにinterface付き。

program test
    implicit none

    interface
        function inv_lapack(A)
            real(8), intent(in) :: A(:,:)
            real(8) :: inv_lapack(size(A,2),size(A,1))
        end function inv_lapack

        subroutine show_dmat(A,fmtf)
            real(8), intent(in) :: A(:,:)
            character(len=*), intent(in), optional :: fmtf
        end subroutine show_dmat
    end interface

    real(8) :: A(2,2)
    real(8) :: AI(2,2)
    integer :: i, j

    A(1,1) = 1.d0
    A(1,2) = 2.d0
    A(2,1) = 4.5d0
    A(2,2) = 6.0d0

    write(*,'("A=")')
    call show_dmat(A, "(f9.3)")
    AI = inv_lapack(A)
    write(*,'("AI=")')
    call show_dmat(AI, "(f9.3)")

    stop
end program test

実行すると、まあうまくいったぽい。

$ gfortran -o test test.f90 -llapack
$ ./test
A=
    1.000    2.000 
    4.500    6.000 
AI=
   -2.000    0.667 
    1.500   -0.333 

QZSの軌道

前回作ってみたプログラムで８月１６日のQZSの軌道を描いてみた。
なぜかatanで計算するとyとzがプラマイ逆になってしまう時があったのでatan2を使って衛星軌道は計算した。なんでかな。
どうもc言語は使い慣れていないのがいけないなと思う。こんなことで時間を使ってしまうとは。

で、QZSを「小松」から見上げたら、こんな感じ。

おお、本当に８の字。

暦計算はGPSと全く同じだった。

放送暦の計算（２）

前に放送暦で衛星座標を計算したらRTKLIBのsatpossと100m以上の差が出てしまった。
（ちなみにQZS）

なんでかなと思っていたら、satposを使うと受信機から衛星までの伝搬時間分が補正されてしまうためだった。（前にメモした話のやつ）
単に時刻決め打ちで衛星座標を計算するなら、eph2posを使えば良い。

計算してみると、

------------------------------------------------------------------
2012/08/16 00:00:00 (x, y, z)
(自前)    -23318436.155264, 18586028.182482, -25160301.817461 
(eph2pos) -23318436.155244, 18586028.182504, -25160301.817463 (RTKLIB API) 
------------------------------------------------------------------

一応計算は合っているぽい。

この微妙な差はEkを計算する際のiteration方法で出るのだろうか？？

ちなみにEkの計算は、ただのiterationと、Newton-Raphson法の２種類試してみたのだけど、収束まで1,2回しか変わらなかった。

------------------------------------------------------------------
static double upd_E(const double E, const double M, const double e) {
    /* return M+e*sin(E); */ /* simple */
    return E-(E-M-e*sin(E))/(1-e*cos(E)); /* Newton-Raphson */
}
------------------------------------------------------------------

放送暦の計算

お盆が終わってから週末がなんだかんだと予定があり、何もできなかった。

今日は早めに帰ったし、久々にプログラムを書いてみることに。
ということでQZSSの放送暦から軌道を計算してみた。

とはいえ、IS-QZSSを読むのが面倒なのでIS-GPSの計算をやってみたのだけども。
(IS-GPS-200E, 8 June 2010 p.92-96)

で、自前で作ったプログラムと、RTKLIBで計算したQZSの座標値を比べると・・・

------------------------------------------------------------------
2012/08/16 00:00:00 (x, y, z)
(自前)    -23318436.155264, 18586028.182482, -25160301.817461 

(satpos) -23318319.272650, 18586114.815274, -25160339.595141 (RTKLIB API) 
------------------------------------------------------------------

なんかずいぶん違うよな・・・。100m超。

相対論効果を考慮してないのだけど、そのせいなのか、はたまた計算ミスなのか。
それとも計算式がGPSとQZSで違うのか。
ただ、酒を飲みながらプログラムを書いたので計算ミスの可能性が非常に高い。
いや、相対論効果って、ここで効いてくるのか？？どんどん分からなくなってきたぞ。
お酒怖い。

しかし、相対性理論が普通の生活で使われていて、こんなに影響してたら、すごいよなあ、と思う。一度まじめに相対性理論を勉強してみたい。

さて、衛星軌道の計算が大体合っていそうなら、あとは最小二乗で位置を計算すれば単独測位プログラムの完成だ。

電子基準点 小松(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コマンドも派生版がたくさんあるようなので、環境毎に対応するしかないのかな。今の美しいスクリプトを汚くするのも嫌だし。

RTKLIBのrtkpostをLINUXで動かす

RTKLIBってwindows版がメインなのか、残念・・・
LINUX版はコマンドラインのappしかなくてマウスでぐりぐりできないよ・・・

て思う人も結構いそうだけど、何気にwineで動くんだよね。

$ wine rtkpost.exe

rtklibを使ってRINEX L1, C1, L2, P2を表示させる

rtklibというGNSS解析ソフトウェアがある。
    http://www.rtklib.com/

なんとフリーで、オープンソース。
ざっと中身を見ると、APIもかなり充実していて、インターフェースが非常にシンプル。単独測位くらいまでならほとんど自前でコーティングせず作れてしまいそう。すごい。

とりあえずRINEXファイルを表示だけするものを作ってみたので使い方をメモ。

rtklibを使って少しずつGNSS解析の中身を理解していきたい。
(*C1, P2と書いたけど、rinex.cの中で優先順位が決まっているので、きちんと固定するには setrnxcodepriを呼ばなければいけない）

 /*
 * printobs.c
 *
 *  Created   : 2012/07/28
 *      Author: sako
 */
#include "stdio.h"
#include "stdlib.h"
#include "rtklib.h"

/* copied from postpos.c in rtklib v2.4.1 */
static int nextobsf(const obs_t *obs, int *i, int rcv) {
    double tt;
    int n;
    
    for (;*in;(*i)++) if (obs->data[*i].rcv==rcv) break;
    for (n=0;*i+nn;n++) {
        tt=timediff(obs->data[*i+n].time,obs->data[*i].time);
        if (obs->data[*i+n].rcv!=rcv||tt>DTTOL) break;
    }
    return n;
}

int main(int argc, char *argv[]) {
    char *rinex_o;
    int i, j, nu;
    char tmstr[20];

    obs_t obs={0};
    nav_t nav={0};
    sta_t sta[1];

    /* parse argument */
    if (argc > 1) {
     rinex_o = (char *) malloc( strlen(argv[1]) + 1);
     strcpy(rinex_o, argv[1]); 
    }
    else {
        printf("no input file\n");
        return -1;
    }

    /* read rinex obsfile */
    if(readrnx(rinex_o, 0, &obs, &nav, &sta[0]) <= 0) {
        printf("Can't read rinex obs file.\n");
    }

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

        for(j=i;j<i+nu;j++) {
            time2str(obs.data[j].time, tmstr, 0);
            printf("%2d: %s %3d %f %f %f %f\n", j-i, tmstr, obs.data[j].sat,
                obs.data[j].L[0], obs.data[j].P[0],
                obs.data[j].L[1], obs.data[j].P[1]);
        }
        i += nu;
    }

    return 0;
}

教科書を注文

GNSS - Global Navigation Satellite Systems: GPS, GLONASS, Galileo and More [ペーパーバック]

Bernhard Hofmann-Wellenhof (著), Herbert Lichtenegger (著), Elmar Wasle (著)

を注文した。

去年の秋頃から始めたGNSS解析の勉強だけど、実は教科書らしい教科書は持っていなかったので、不便で仕方がない。

8,198円・・・。もう注文ボタンは押した。


やっぱりデータをプロットして直接目で見てみないとイメージがわかないので、少し解析プログラムを組んでみたい。


widelane+narrowlaneくらいを勉強がてら自分で計算してみたいけど、基準衛星の選択とかが面倒くさそうな気がする。座標変換も面倒くさそうな気がする。

ああ、こんなこと考えながらwidelaneてどんな式だっけか、と今日思ったのがこの教科書を買う動機だったのを思い出した。


と、理想はあるけど、今のところ書いたプログラムはまだRINEXの読み込み、編集、MP線形結合しか計算できない。全然役に立たないな。

そういえば、teqcで国土地理院から出たGPS+GLONASS+QZSのデータを結合しようとしたら、エラー。

teqc: failure to read '     2.12           OBSERVATION DATA    M (MIXED)           RINEX VERSION / TYPE' on line 1 of './obsdata/05481951.12o'
    (unaccepted RINEX version or non-RINEX file; must be RINEX Version <= 2.11) ... exiting

ええ！teqcってRINEX ver.2.12に対応してないのか！？
じゃあ作ったプログラムはファイルの結合と、ヘッダにプログラム名を入れて自己満にひたるために使うか。

GNSS IFB

ブログをはじめることにした。

GNSSの勉強やプログラミングのメモなどを書いていこうと思っている。

今日読んだのはこれ。

Jean-Marie Sleewaegen, Andrew Simsky, Wim De Wilde, Frank Boon and Tom Willems

Demystifying GLONASS Inter-Frequency Carrier Phase Biases

 http://www.insidegnss.com/node/3073


これまでGLONASSを別の種類の受信機で二重位相差を取ると発生する周波数依存のバイアスは、回路の群遅延だと思っていた。

でもこの論文によるとその影響はせいぜい0.1mm程度で、数cmにも達するこのバイアスは受信機ファームウェア内の信号処理に起因するとのこと。

よって、温度にも依存しないし、一定値となることが期待される。
これが本当ならバイアスの傾きを定数として求めておくことは十分可能だと思われる。

しかも、周波数に一つ、ではなく、受信機に一つだけ求めておけば良い。
(距離に直すと同じ値のバイアスになるため)

・・・しかし、これは、本当に本当なんだろうか。

テスト

テストです

ブログをはじめることにした。

GNSSの勉強や、プログラミングのメモなどを書いていこうかと思っている。