私はGNSSを10年来、仕組みをよく知らずに使ってきました。しかしUAV写真測量でも、GNSSによる干渉測位は標定点の測量、RTKドローンによる撮影位置の測位に広く使われていて、精度の要であるともいえます。UAV写真測量の高精度化・効率化を追求するなら、GNSSの基礎知識も必要に思われます。
そこでついに最近、何事にも重い腰を上げ、高須先生のオープンソースソフトウェアRTKLIBを手掛かりに、GNSS測位の勉強を始めました。
写真測量におけるTriangulationを既知点（カメラ）からの「向き」に基づく測量と呼ぶなら、GNSS測量は既知点（衛星）からの「距離」に基づく測量であり、
両者には似ているところもありそうです。
この「GNSSの勉強メモ」シリーズの記事は、勉強でわかった／わからないことのメモ書きで、随時、加筆・更新していくような、勉強ノートにするつもりです。私自身にわかりやすいように書いているだけですが、もしかすると、同じところで躓いた方の参考になるかもしれない、または逆に詳しい方に何か教えてもらえるかもしれないと考え、ここに書くことにします。
初回はRTK-GPSの仕組みに関するメモですが、RTK-GNSSと書かないのは、まずは米国のGPSだけを使ったシンプルな場合の原理を勉強しているためです。
（色の凡例： 緑字は自信のない予想（未確認） 赤字は疑問点（未解決））
- 高須 (2007)によると、RTK-GPS測位の基礎となる方程式は式(A.4.4)だ。基線が十分に短いという仮定が許される限り、本質的にはこれを多数の衛星の組合せについて連立し、未知数、つまり移動局の受信機の位置と整数バイアスを求めるだけのようだ。
- ところが、式(A.4.4)の第1式に含まれる「二重差整数バイアス」ともよぶべきN_ab_urは未知であり、衛星の組ごとに変わる。従ってこの式を何組の衛星について立てたところで、方程式の数は未知数の数に追いつけないため、解は不定となる（※1）。
- 実は式(A.4.4)は2式からなり、1つめは搬送波位相の観測値に関する方程式だが、2つめは擬似距離の観測値に関する方程式である。一見、この2つめがあるおかげで、方程式の数は未知数の数に勝つことができる。例えばm個の衛星と1周波の受信機を使うと、独立な方程式は2×(m - 1)本、未知数は移動局受信機の位置について3、二重差整数バイアスについてm-1だから計m+2個であって、m≧4のとき、方程式の数≧未知数となる。m>4のときは方程式が余るので、最小二乗法などで統計学的に合理的な解を求められる。
- しかし以上の解釈によれば、干渉測位なのに、擬似距離が主力の一翼を担っているように見えてしまう。それでは精度が出ないはずなので、おそらく実際には、式(A.4.4)の第2式は、第1式に比べて小さな重みが設定されているはずだ。実際、いま勉強しているRTKPOSTにも"Code/Carrier‐Phase Error Rate"という設定項目があり、その2乗の逆数くらいが第2式の重みとして使われているのではないかと予想する。また、第2式は受信機と衛星との間のおよその距離を示してくれるので、後述のカルマンフィルタによる解法において、解を狭い範囲に効率よく導く役割があるのではないかと思う。
- 擬似距離の方程式にあまり頼れないなら、方程式の不足を解消する本命の策は何だろう。それは、式(A.4.4)を複数時刻（複数エポック）について連立することのようだ。A.4.2.3のようにカルマンフィルタを使うということは、結局、そういうことのはずだ。私はカルマンフィルタに馴染みはないが、こちらのページに書かれているように、これはモデル（方程式）に基づいて何かを推定する方法であり、ここでは最小二乗法と同じく、連立方程式の解法に近いものと見なせるからだ。解法として最小二乗法を示している他の資料でも、複数時刻に関する搬送波位相の方程式の連立を用いている分かりやすい資料を複数見つけた（浪江 (2001)や、RTKの話ではないが舞鶴高専の四蔵先生の講義資料）。
- 複数時刻について連立するとき、未知数のうちのどれを一定として設定するかが重要だ。すべての未知数が時刻間で独立であるとみなしてしまったら、連立するうまみがないし、方程式の不足も解消できない。まず、移動局の受信機の位置について、スタティック測位ならば、移動局（※2）の受信機の位置は一定（全方程式で共通）とみなせるが、RTKやPPKの場合はそうはいかない。高須 (2007)のA.4.2.3のp.4左上にあるカルマンフィルタの時間更新則の説明によると、そもそも時間更新の時点で全く独立になるとして扱われている：時間更新則において、受信機の位置はその時刻の単独測位値に置き換えてしまい、整数バイアスだけを引き継いでいる（※3）。次に整値バイアス（※4）については、これはサイクルスリップがない限り、観測（位相積算）を続けていれば変わらないので、RTKでも前文の通り引き継ぐことができる（p.4左上のF_kの式参照 ※5）。
- 以上のような連立方程式を最小二乗法などで解いたとき、「整数」バイアスとは名ばかりでその解は実数値であり、このときの測位解（移動局の受信機の位置の解）はFLOAT解と呼ばれる。この時点では、整数であるという拘束条件をつけていない（※6）のだから、当たり前の結果だ。整数バイアスが整数であるという拘束も利用した方が測位解の精度が向上するので、その後、実数値の解に「近い」整数値を整数バイアスとして選び、晴れて整数となった整数バイアスに基づく測位解：FIX解を得る。
- この整数バイアスの「決定」作業が、RTK, PPK, スタティック測位などにおけるもう1つの山場のようだ。二重差整数バイアスは衛星の組ごとにあるので、例えばm個の衛星と1周波の受信機を使う場合にはm-1個ある。それぞれの衛星の組について独立に、四捨五入によって一番近い整数を選べばよいなら単純だが、それだと元々の連立方程式の解として適さない整数バイアスの組になる可能性がある：最小二乗法の場合、残差平方和の大きい解になる恐れがある。対策として浪江 (2001)のp.10には、可能性のある整数バイアスの組を絞った上で、残差平方和を最小にする組み合わせを採用する決定方法の記述がある。浪江 (2001)がその前に説明している解法では、各時刻において、最初のエポック(0)との連立によってFLOAT解と対応する整数バイアスを得ているため、時間の経過とともに整数バイアスはある整数付近に収束する傾向を示している（図1.6）。従って、ある程度時間が経てば、確かにあり得る整数の組を少数に絞ることができそうだ。
- 一方、高須 (2007)は式(A.4.12)で、浪江 (2001)とはおそらく別の、整数バイアス決定方法を述べている。この式については、直上の説明中の表現「実数推定誤差共分散により定義された距離」の解釈が難しいが、RTKLIB 2.4.2のユーザーマニュアルの式(E.7.17)と同じ意味だと思われる。式(E.7.17)によれば、選ぶ整数の組は「（カルマンフィルタで実数値として）推定済みの整数バイアスとの差の重み付き二乗和が最小になる」ような組である。つまり、単純に言えば推定済みの実数値と近いものを選んでいる。ただし、推定したときに得た推定値の共分散行列Q_Nの逆行列で重みづけしているので、推定値の分散が大きい衛星の組の整数バイアスは軽視される（推定済みの実数値と選ばれる整数の差が大きくても許容される）し、推定値間の相関も考慮されている（例えば、推定値の相関の強い衛星の組がある場合、相関がない場合と比べると両方とも軽視される）。これを浪江 (2001)と比べると、「推定値の実数値に近い整数を選ぶ」ことがより徹底される一方で、「元の連立方程式をできるだけ満たす整数の組を選ぶ」という要素がない気がするが、私の理解は合っているだろうか。
※5 RTKLIB 2.4.2のユーザーマニュアルのp.167にあるように、現在のRTKLIBでは、エポックごとに搬送波の一重差整数バイアスをリセットする設定：Integer Ambiguity Resolution = Instantaneous もある。しかしこうしてしまうと、擬似距離の観測値に関する式(A.4.4)の第2式に頼らないと必要な方程式数を確保できなくなるので、精度は落ちるのは必定に思われる。
※6 なぜ最初から整数条件を付けないのか。私の予想では、(A.4.8)のような未知数による微分ができなくなるからだろう。最小二乗法を使う場合でも、残差平方和を未知数で微分できないと、「残差平方和の未知数による微分＝0」の点として解を求めるアプローチが使えず、広大な空間の探索に頼ることになってしまいそうだ。