logo 質問と回答です
ホームコンパイルの方法BCCの使い方質問と回答記事の正誤表履歴電子研ホーム


【4】高須知二様(技術コンサルタント)から頂いたご意見に関する捕捉
ご意見の趣旨 (1)GPS衛星位置を求めるため,信号送信時のECEF座標系から信号受信時のECEF座標系に変換する際, OMEGADOTE*p_range/(C)だけZ軸回りに回転を行っています. p_range/(C)は信号の伝搬遅延を表していると思いますが,擬似距離p_rangeには,衛星時計及び受信機の時計誤差分が含まれます ので正確な信号伝搬遅延分にならず,これらの誤差は無視できません.
(2)「質問と回答」にあるomegadotへ乗じる時刻の件ですが,ICDを忠実に解釈するという意味では,元の式が正しいと思います. 放送歴の精度を考えると実効上差はないのですが,この修正でノバテルと一致性が高いとするとノバテルの計算には問題があるかもしれません.

補足の解説 貴重なご意見頂きまして大変ありがとうございました.(1)につきましては,解説記事中での説明を省略しておりました. ご指摘はごもっともと思いますので,以下に捕捉させて頂きます. 衛星時計の誤差分につきましては,プログラム中で,af0,af1,af2を使った二次式で擬似距離(psr)を補正しています. 具体的には,calc_sat_position.cの中で関数eph2xyz()をよぶ前の部分になります. 次に受信機時計誤差分ですが,公開しているプログラムでは補正しておりません. ご指摘のとおり,通常の受信機(Trimble,Ashtech等)では, 擬似距離に対する受信機時計誤差の補正がある一定時間間隔で行われますので,受信機時計誤差分が無視できず, 光路差方程式を解いて補正する必要があると思われます. ノバテルの受信機では,擬似距離に対する受信機時計誤差の補正が毎エポック行われています.このため,ノバテル受信機に限っては,光路差方程式を解く補正は必要ないと考えました. また厳密にICDに従えば,送信時刻を求めるために,相対論項とT_GDの補正が必要ですが,小さいので無視しています.
(2)につきましては,再度考え直し,次段落の質問の回答を少し訂正しました. 確かに放送歴の精度以下ですので,大きな問題はないと思いますが,少し疑問が残ります.
【3】楊海龍様から頂いたご質問(ご助言)とその回答(解説)
Q.(その2)でGPS衛星位置の計算結果とノバテル受信機の結果の差がX,Y軸で2cmになる問題(Z軸は1mm)ですが,いろいろと考えてみました.この結果,プログラムを以下のように書き直し,1時間データで確かめたところ,X,Y軸についても1mm以下となりました.これは,計算式の一部を送信時刻から受信時刻に変更したためです.この理解は正しいでしょうか?また,何故このようになるのでしょうか?

(変更前)
omegak = omega0 + (omegadot - OMEGADOTE) * tk - OMEGADOTE * (toe + p_range/(C))

(変更後)
omegak = omega0 + omegadot * (tk + p_range/(C)) - OMEGADOTE * (toe + tk + p_range/(C))

A. すばらしいご助言を頂きまして,ありがとうございました.私の方でもこの問題についてもう一度考え直してみました.私の解釈の誤りがあるのかもしれませんが, 以下に変更後の式を採用する理由を考えてみました. まず,(その2)の最後で述べましたように,「衛星位置は,信号の送信時刻で計算され,かつ信号が受信される時刻のECEF座標での位置」でなければいけません. そこで私の計算では,昇交点赤経(omegak)を求める部分(プログラムcalc_sat_position.c:関数eph2xyz())で, 信号の伝搬時間(p_range/(C))だけtk(送信時刻)を進め, 地球の自転レート(OMEGADOTE)に乗じることで,地球の自転を戻し,受信時刻のECEF座標での衛星位置を求めておりました(変更前). (変更後)の式では,omegadotに乗じる時刻も送信時刻tkではなく,受信時刻にしています. これは,昇交点赤経(omegak)が,週始め(tow=0)の昇交点赤経(omega0)と,omegadot(昇交点赤経の補正レート)と toeからの経過時間の積で求める補正項の和により計算されるため, この補正項についても,送信時刻tkではなく,受信時刻を使うという解釈と思われます. (第3項はomega0が週始めだから) 変更後の式を用いて24時間データで確認したところ,ノバテル受信機の計算結果との差は,Z軸だけでなく,X,Y軸についても最大1mmであることが確認できました.

【2】上島一彦様から頂いたご質問#2とその回答
Q.最小2乗法の重みに関してですが,本文では「重みは,誤差分散に反比例するように付ける」と記述があります.しかしプログラムリストでは,分散(単位はm2)でなく,rng_std(単位はmなので標準偏差)に反比例するように付けています.これは何故でしょうか?

A.本文中6.で簡単に説明しているのですが,わかりにくいかもしません。このプログラムでは,最小2乗法にハウスホルダー法という手法を用い,観測方程式から直接,最小2乗解を求めています.これは,5.で説明している,式(11)とは異なる解法で,式(5)または(7)に重みづけしたものを直接解いています.式(7)の重み付けについては,6.の最後の辺りに記述があるのですが,式(5)のΔRをW1/2ΔRに,AをW1/2 Aに置き換えます.これは,式(5)を重み付けで書き直したとき,

  W1/2ΔR = W1/2 A X

となることを意味します.したがって,プログラム中のmake_weight()で求めているのは,WではなくW1/2(対角要素に1/σが入った行列)で,lls()に渡している行列(matrix_aw)は上式の右辺のW1/2 Aであり,行列(delw)は左辺のW1/2ΔRです.

以上から,make_weight()で計算するのは,W1/2の要素で1/rng_stdとなります.この重み付きの観測方程式については,文献[3.2]p.3-37に詳しい説明があります.
【1】上島一彦様から頂いたご質問#1とその回答
Q.数値計算ハンドブックの Pro.8.1.1 を調べてみました。ところが、Pro.8.1.2 の行番号10で、暗礁に乗り上げてしまいました。
   10 A(I,J)=CC-DFLOAT(ICC)
突然現れた変数CCとICCの意味が不明です。

A.私の方でハウスホルダー変換をチェックしたプログラム(Pro.8.1.1)を調べましたところ,明らかに本の記述とは異なり,ご質問の行の前に,
   ic = ic + 1;
   cc = (double)ic * c;
   icc = cc;
が追加されています.たぶんこれは,本の誤植だったと思います.どうやって調べたかは,ずいぶん前のことでもあり,記憶がはっきりしないのですが,この本の旧版に正しいプログラムが掲載されていたか,または出典:REFERENCE...1965(Pro.8.1.1の最初のコメント行にある)を調べたのかもしれません.それで,このプログラムを使って,サンプルの実効結果(表8.1.1)が得られることは確認できています.





enri_print
For comments or suggestion, please contact me.