空中測量研究室の技術ノート【2冊目】

山口大学の1研究室による研究メモです。UAV写真測量, ドローン測量, フォトグラメトリ, SfM/MVSなどと呼ばれる技術の情報があります。

【GNSS勉強メモ】SNRとCNRの違いは何か

GNSSの信号の強度あるいは品質の指標として使われる、SNR (Signal-to-Noise Ratio; S/N)CNR (Carrier-to-Noise Ratio; C/No)の違いは何だろうか。

私が少し調べた限りでは、

  1. 人により、同じ意味で使われることもあるようだし、定義が異なることもあるようだ。
  2. Richardson et al.(2016)によると、SNRとCNRはともにコード変調された信号における信号パワーとノイズパワーの比であるが、SNRは受信機内の相関器で計算されるもので、CNR(C/Noと表記されることもある)は(増幅前の)受信アンテナにおける値とのことだ(定義αとする)。信号もノイズも、アンテナと相関器の間で同じくらいの倍率で増幅されるので、CNRとSNRの値はほぼ同じと考えて良いとのこと。
  3. Richardson et al.(2016)に引用されているLangley (1997)によると、C/Noは Carrier-to-Noise density Ratioの略で、ベースバンドの1 Hz帯域幅における搬送波のパワーとノイズパワーの比とされている。そして文脈を見ると、どうやら復調後に評価するようだ。CNRの名前からするとこちらの定義の方が納得がいく(定義βとする)。
  4. IGSによるRINEXフォーマットの定義では、5.7節のTable 12などを見る限り、S/NとCNRは同義に使われている。そしてObservation codeがSnaの形式になっているデータは、搬送波の観測がある限り、搬送波に関するものつまり上記βに対応するもののようだ。例えば下図の例では、データ部の4列目がこれに該当し、最初のエポックにおけるG21のCNRは45だ(単位はおそらくdBHz)。なお、RTKLIBのRTKPlotではRINEXファイルのこの列をそのままSNRとして表示しているようだ(RTKPlotとテキストエディタで同じ.obsファイルを開いて確認)。受信機を制御するアプリDrogger GPSでリアルタイムで見られるCNRや、NMEAログに記録できるCNRがこれ(RINEXファイルで見られるCNR)と同じかどうかは未確認
  5. RINEXファイル (ver.3.02)の冒頭部分の例
  6. 定義βにおけるCNRは、マルチパスにどのように影響されるのだろうか。直達波に、波形が崩れた反射波が重なっただけなら、ノイズが増えたことになってCNRは低下しそうだ。一方、直達波に、波形を保った反射波が重なってしまったら、位相の関係によって搬送波の振幅が見かけ上大きくなったり、小さくなったりするはずだ。従って一概に、CNRが高ければマルチパスの影響が小さいとも言えないように思われる。以前の記事でも引用させてもらった鈴木 (2016)によれば、「マルチパスの影響によってその値が落ち込んだり、また逆に増幅されたりなど値が大きく変動する。しかし、マルチパスが存在しないときには、C/No 比は衛星の仰角による。 」とある。またこちらの記事によると、マルチパスは誤差のうち比較的高周波の変動成分の原因とみなされるようだ。つまり、例えば各衛星のCNRが仰角の滑らかな関数のように推移し、小刻みで大きな変動がなければ、マルチパスの影響は少なそうと言えそうだ。

 

【GNSS勉強メモ】水面によるマルチパス関連

UAV写真測量はGNSSにべったり頼るものだが、私のGNSSの勉強は遅々としてあまり進んでいない。次に勉強する時間が手に入ったときのために、現状の疑問点・課題をメモしておく。赤字は今後の課題となるところ。

 

  • 水面は、マルチパスの原因の1つとして目されているらしい。広い水面の見える場所にアンテナを置く場合には注意が必要のようだ。
  • マルチパスの影響には、擬似距離の計算に対する影響と、搬送波の位相測定に対する影響があるらしい。
  • 擬似距離の計算において、マルチパスは「受信機内のコードのレプリカと、受信電波に乗ったコードの相関をとって時間差を測る処理(相関がいちばん強くなる時間差を求める処理)」を邪魔する働きをする:
    https://www.denshi.e.kaiyodai.ac.jp/wp-content/uploads/pdf/investigation/030727.pdf
    https://www.denshi.e.kaiyodai.ac.jp/wp-content/uploads/pdf/investigation/20040523.pdf
    が、反射波と直達波の伝搬距離の違いが大きければ、相関のピークとなる時間差への影響が小さくなるので測位解に影響しにくい。
  • 水面反射波と直達波を比べたときの伝搬距離の増分を、平面波を仮定して幾何光学的に考えてみると、下図のようになった。数値例のとおり、アンテナの水面からの高さが数 m程度と低くては、伝搬距離の差もコードの波長に比べて小さなものになる。また同じアンテナ高さなら、入射角が大きい反射波ほど、つまり水平線近くから飛んでくる反射波ほど、伝搬距離の差が生じにくく、コード測位にとっては厄介そうだ。

  • 一方で搬送波については、反射波が重なることで位相が変わるので、そもそも伝搬距離の違いが大きくても搬送波位相測定への影響が小さいとは言えないように思う(空気中での減衰を無視すれば)。水面の各部からの反射波のうち、無害な反射波とは、上図の2 h cosθが大きい反射波ではなく、それが波長で割り切れる反射波か、水面での反射率が小さかった反射波だろう。
  • 水面を平面とみなしたとき、電波は水面で鏡面反射することになるが、その反射係数や反射率は、私がよく使ってきた可視光と同じくフレネル (Fresnel)の式に従うようだ

屈折率1の空気から屈折率1.33の水に可視光が入射するときの反射係数(フレネルの式)

屈折率1の空気から屈折率1.33の水に可視光が入射するときの反射係数(フレネルの式)
  • ただし、Fresnelの式に登場する水の屈折率は周波数に依存し、光(約1.33)と電波では異なり、Fresnelの式を使いたければGNSSで使われる電波の周波数における屈折率を調べるか、屈折率=√(比誘電率×比透磁率) で計算する必要がある。
  • こちらの論文には、GPSのL1波の周波数における0℃の淡水の比誘電率(Table 1)や、屈折率(Table 2)が示されている。屈折率は9.24のようだ。それをFresnelの式に代入して描ける反射係数VS入射角のグラフは下図のようになる。この論文のFig. 2は、p波、s波ではなく、co-polarizedとcross-polarizedの反射係数を示しているので、下図とは見た目が異なる。ちなみにしっかり調べていないが、比透磁率は可視光のみならずGNSSの周波数においてもほぼ1とみなされているようだ。

屈折率1の空気から屈折率9.24の水にL1波が入射するときのフレネル反射率

屈折率1の空気から屈折率9.24の水にL1波が入射するときの反射係数(フレネルの式)
  • なおこちらの論文の式(4)は、屈折率の代わりに比誘電率が√なく使われていて、Fresnelの式と異なるように見える。こちらの別の論文でもそうだ。なぜだろうか、要勉強。
  • こちらの論文のFig. 2では、入射角が増加してBrewster角に達したときにco-polarizedの反射係数がcross-polarizedの反射係数を上回るが、この時点で水面に入射する右旋円偏波(RHCP; GNSS衛星では右旋円偏波を使用)が左旋円偏波 (LHCP)に切り替わるらしい。反射波のp波とs波の振幅が違うのに、だ円でなく円偏波になるとはどういうことだろうか。円偏波まわりの勉強も足りず、このあたりの意味がまだ分かっていないが、こちらの別の論文のpp.39 - 40にも同じようなことが書いてあって参考になる。式(4.11)(4.12)を使うと、L1波の反射係数は次のようになる。ただしこちらの論文Fig.2と合っておらず、何か勘違いしている可能性がある。

右旋円偏波が入射した場合の反射係数

屈折率1の空気から屈折率9.24の水にL1波(右旋円偏波)が入射するときの反射係数(フレネルの式)。Brewster角以下では式(4.11)、以上では式(4.12)を使用。
  • GNSS用のアンテナでは、RHCPのゲインを高く、LHCPのゲインを低くするように作られているので、Brewster角よりも大きい入射角で水面の入射した波については、反射波が左旋円偏波になるため、心配が少ないということになるようだ。
  • こちらの別の論文の表4.1によると、海水 (Sea Water)と淡水 (Fresh Water)では比誘電率が大きく異なる。淡水が80なのに対し海水が20となっている。これがL1波の周波数における値だとすると、海水の屈折率はおよそ√20=4.5くらいになるだろうか。気にしている水面が海の場合は、上のグラフは該当しないので注意。
  • マルチパスに強いGNSSアンテナというのは、(グランドプレーン、チョークリングの使用も含めて)水平面下からの電波のゲインを抑えられるアンテナや、「軸比」というものが1に近く左旋円偏波の影響を受けにくいアンテナのようだ。
  • 軸比 (Axial Ratio)の定義を調べるのに苦労した。その理由として、アンテナの「利得(ゲイン)」と似て、検索すると送信アンテナと受信アンテナの話が両方ヒットすることもあるが、定義が複数あるようだ。こちらの資料の式(2・58)と、こちらの資料のp.33では、符号が逆になる式が示されている。おそらくGNSSの文脈では後者か。よく「長軸と短軸の比」と書かれているが、何の楕円偏波の話をしているのかわからない。はっきり定義がわからないと、仕様として書かれた軸比の最大値(今使っているアンテナでは3 dB)から、水面反射波のゲインのの最大値を計算することなどはできない。
  • RTKLIBのRTKPlotでもマルチパス [m]の表示ができるが、私にはその定義はわかっていない。RTKLIBとは別に、マルチパスの解析のためのオープンソースのソフトがあって、そのページにマルチパスの推定式が載っていた。これを手掛かりに、このような足し算で何を推定しているのか、物理的な意味などを調べていきたい。

 

Metashapeでのカメラパラメータ・タイポイントの分散の調べ方

1. はじめに:最小二乗法による推定値のバラつきとは?

何かのパラメータを最小二乗法で推定するときは決まって、モデル式に(多くの場合、モデル式に含まれる観測値に)何かの誤差が含まれていると仮定している。もし誤差を認めないなら、つまり誤差項のないモデル式がすべて厳密に成り立つなら、必要最小限の数の式を適当に選んで、連立方程式として解けば済むからだ。

 

そしてその誤差に特別な仮定を置けば、求めるパラメータの推定値とともに、その推定のバラつきの指標(分散あるいは標準偏差)を得ることができる。一番簡単な例として、n組のy, xのデータから単回帰直線を推定する状況(単回帰分析)では、

y = ax + b + ε

というモデル式を考え、最小二乗法でパラメータa, b(直線の傾き・切片)を推定する。ここで、誤差εは、yの観測値に含まれる誤差を表している。そしてこのn個の観測誤差が、互いに独立で、同じ1つの正規分布に従うと仮定すれば、a, bの推定の分散を計算することができる(この単純な状況での分散は、比較的簡単な式で表されるので、調べてみよう;また以上の説明がピンとこないなら、Excel擬似乱数を発生させる関数を使い、自分でx, yのデータを作ってa, bを推定してみよう)。

 

この分散とは何だろうか。推定値は1つしかないのに、なぜバラつきがあると考えられるのだろうか。答えを書くと、この分散は、観測誤差に関する上記の仮定のもと、「観測データの生成とそれに基づく最小二乗法によるパラメータa, bの推定」を無限回行ったときに、結果として得られるa, bの推定値のバラつきを表している。誤差項の確率分布について仮定を置くことで、手元にあるデータ自体を「いくらでも繰り返せる試行の中の、ある試行で生成されたデータ」とみなすことができ、同様の試行を無数に考えることができるわけだ。

 

話変わって、セルフキャリブレーション付きのSfMあるいはバンドル調整では、内部パラメータ、外部パラメータ、タイポイントの3次元座標を同時に推定する。推定するパラメータがa, bの2つしかなかった単回帰分析の場合と比べて、モデル式も複雑だ。それでも、タイポイント(元は特徴点)の画素座標(u, v座標)の観測に誤差を認め、それが互いに独立で、1つの正規分布に従うと仮定すれば、単回帰分析の場合と同様に、それぞれのパラメータの推定値について、その推定の分散を推定することができる。最近のMetashapeでは、それらをかなり見ることができる。厳密な話をすると、Metashapeが表示する分散が、誤差について上記の最も単純な仮定を置いた場合の分散なのかどうかは未確認であるが、この記事では、それらの分散のチェック方法をまとめてみる。

 

2. 内部パラメータの分散の確認方法

内部パラメータの分散は、その平方根である標準偏差の数値として、"Camera Calibration"(カメラキャリブレーション)ダイアログから表示できる"Distortion Plot"(歪曲プロット)の"Correlation"(相関関係)タブの、"Error"(誤差)という名前の列で確認できる。

Metashapeの相関関係ダイアログ

Metashapeによる内部パラメータの相関係数行列・標準偏差の表示

列名の表示が"Standard Deviation"(標準偏差)などでなく、"Error"(誤差)となってしまっているため、本当に標準偏差なのか、この画面を見るだけでは確信が持てないが、ユーザーマニュアル(ver. 1.8.4のもので確認)には、これが標準偏差であることが明記されている。

またこのタブには、他の記事でも触れているが、内部パラメータ間の相関係数が行列(相関係数行列)として表示されている。相関係数行列と上記の標準偏差を使えば分散共分散行列が計算できる(2変数の相関係数に両変数の標準偏差を掛ければ共分散になる)ので、内部パラメータについては一応、Metashapeでは分散共分散行列が取得できるとも言える。

 

3. 外部パラメータの分散の確認方法

外部パラメータの分散は、標準偏差の数値として、"Reference pane"(座標データペイン)で確認できる。具体的には、"View Errors"(分散を表示)というボタンを押すことで、カメラの投影中心の座標を表す3列と、回転角を表す3列のそれぞれについて、列名の末尾に"sigma"(シグマ)が付いた列が登場する。"sigma"は標準偏差を表すので、「分散を表示」という訳のボタン名ながら、実際に表示されるのは標準偏差のようである。

外部パラメータの標準偏差を表示

外部パラメータの標準偏差を表示

 

4. タイポイントの分散の確認方法

タイポイントの分散は、アラインメント後のオンデマンドのバンドル調整を行うダイアログ"Optimize Camera Alignment"(カメラアラインメントを最適化)で、"Estimate tie point covaiance"(タイポイントの分散を推定する)にチェックを入れることで、計算できる。

タイポイントの分散を推定

"Estimate tie point covaiance"にチェックを入れてバンドル調整

推定した分散を表示するには、Model(モデル)メニュー -> View Mode(ビューモード) -> Point Cloud Covariance(点群の分散)を選択すればよい。Modelビューにベクトルプロットとして表示される。これはユーザーマニュアルによると、それぞれのタイポイントについて、誤差楕円体 (error ellipsoid) の長軸の方向と長半径(長軸の長さの半分)を示したものである。色も付けられていて、分散が大きいタイポイントは赤く色づけられている。

 

Modelビューで、タイポイントの分散を表示

Modelビューで、タイポイントの分散を表示

 

これは何かというと、上記のような無限回の「試行」をしたとき、タイポイントが最も動きやすい方向とその動きやすさを示したものである。タイポイントの座標にはX, Y, Zの3成分あるので、X, Y, Zそれぞれに分散があるし、X-Y, X-Z, Y-Zそれぞれの間に共分散がある。それらの兼ね合いで、最も動きやすい(バラつきの大きい;不安定な)方向と動きやすさ(バラつきの大きさ;不安定さ)が決まる。このベクトルプロットは、それらの視覚的に示したものである。詳しく学ぶには、まずは2次元の「誤差楕円」について学ぶのがよいと思われる。

 

タイポイントは3次元復元したい被写体上にあるわけだから、タイポイントの分散は、検証点誤差など対空標識が必要な指標を除けば、3Dの成果物の不安定さに関するもっとも直接的な指標になるだろう。例えば、SfMの後にMVSをして密な点群を作るなら、その点群の分散は、SfMによって生じるこのタイポイントの分散に、MVS自体で生じる分散が乗ったようなものになるはずだ。つまりタイポイントの分散は、被写体の3Dモデルの各点の位置のバラつきに直結する。このベクトルが長い領域では、検証点誤差やMVSで作る点群の位置の誤差が、小さくなることは期待できない。このベクトルが対象領域全体にわたって長ければ、それは撮影計画が悪く、内部パラメータの一部の分散が大きい(不定に近い)可能性がある。例えばfが不定に近く、fのみが不安定な状況があったとしたら、鉛直方向の長いベクトルが観察されると予想される。

 

なお、「タイポイントの分散」があれば検証点誤差の評価は要らないかと言えば、次のような理由から、そんなことはない。

  • 分散は誤差のバラつきの指標であって、誤差の平均(バイアス)の評価を含んでいない。非線形モデルについて、最小二乗法による最尤推定量が不偏推定量と言えない限りは、バイアスの評価も必要なので、検証点誤差も必要だ。
  • 分散は、無限回「試行」したときのバラつきを示すものであって、たまたまその1試行である現実が、どのくらいの誤差をもっていたかを示してくれるものではない。タイポイントの分散が小さくとも、運が非常に悪ければ、大きな検証点誤差を生じることもあり得る。

 

Metashapeを使ったCGシミュレーションでCritical Configurationを調べる

本記事では、前記事:バンドル調整でパラメータの不定性を確認するにはどうするか を踏まえ、前記事で第2の検証方法として挙げた「2. 最適化を実施したときの目的関数値が、θiを誤った値に固定した場合とθiを正しい値に固定した場合で、有意に変わらないかを調べる。」を実演したい。

具体的には、CG空間上で、完全な水平面を、理想的なピンホールカメラを鉛直下向きに向けて撮影した画像セットをMetashapeで解析して、この「平行光軸撮影で、光軸に直交する平面状の被写体を撮ること」が、内部パラメータの一部にとってCritical Configurationであるかどうか(不定になるかどうか)を、調べてみる。不定性を調べる内部パラメータとしては、次の2つをとりあげる。

①準備

使用するCG画像はこんな感じだ。対空標識も写っているが、平面上にあり、飛び出してはいない。6年ほど前のCGシミュレーションに使った画像を利用している。

CG画像の例

まず、内部パラメータに真値を与え、固定してしまう。理想的なピンホールカメラで撮影(レンダリング)しているので、f以外のパラメータは0だ。

内部パラメータを真値に固定する様子

そしてアラインメントをする。言い遅れたが撮影位置は、このように3×5のグリッド状に配置されている。

アラインメントの結果

アラインメント後に一応、ダメ押しのバンドル調整「カメラアラインメントを最適化」をする。

ダメ押しのバンドル調整

その後、再投影誤差RMSをチェックする。撮影位置・向きの情報も、標定点も使っていないバンドル調整なので、目的関数は再投影誤差RMSだ。最小化されているのは、括弧の外にある無次元の値のはずで、今回は0.0112334だ。これが、fもK1も正しい場合の目的関数値として、ベンチマークになる。

再投影誤差RMSの確認

焦点距離fの不定性の検証

準備ができたので、fの不定性の確認に進もう。焦点距離fを固定したまま、値を誤った値1234に変更する。

fを固定したまま、誤った値に変更

このとき、再投影誤差RMSは大きくなるが、当然だから気にしない。fを固定したまま、バンドル調整を行う。

fを誤った値に固定したまま、バンドル調整

バンドル調整後に、目的関数(再投影誤差RMS)をチェックすると、0.0112335となっていて、fが真値の場合の0.0112334とほとんど変わらない。0.0000001だけ大きいのが気になるが、実はMetashapeでは通常のアラインメント直後の「カメラアラインメント最適化」でも最後の桁が増えることはよくあるので、意味のない差と思われる。さらに、今度はfの値を別の誤った値123, 2345, 3456に固定してバンドル調整すると、いずれの場合もぴったり0.0112334になった。

fの値が誤っていても、最適化後の目的関数値が変わらない。

つまり、fが誤っていても、他のパラメータ(今回の場合、f以外の内部パラメータは全て固定されているので、外部パラメータしかないが)がそれに連動して誤る(fの誤りをカバーしようとする、と言えるかもしれない)ことで、目的関数の増加を防げたということだ。これがまさに、上で言うところのfの不定性を示している。

 

③K1の不定性の検証

次にK1の不定性を調べるため、fに真値2000を入れてバンドル調整をやり直し、②の開始時の状態(①の終了時の状態)に戻す。その後、K1に謝った値-0.05を入れて、

K1を固定したまま、誤った値に

すべての内部パラメータを固定したまま、バンドル調整を実施する。

誤ったK1含め、すべての内部パラメータを固定して、バンドル調整

モデルビューで、点群とカメラの配置双方に、ボウル状の変形が確認できる。K1が誤った値に固定されている条件で目的関数を小さくするため、外部パラメータが連動して誤った(K1の誤りをカバーしようとして誤った)ためだ。

ボウル状変形が発生

しかし興味があるのは、目的関数である再投影誤差RMSだ。今回は0.110368と、K1が真値の場合の値0.0112334より増えている。K1を+0.05にすると0.0798123, -0.1にすると0.542339、-0.01にすると0.0118658と、K1の値によって差は異なるが、真値の場合と比べて一貫して増えている

目的関数が増えている

これは、K1はfと異なって、不定ではないことを示している。つまり今回のconfiguration(撮り方と被写体の組み合わせ)は、fにとってはcriticalであるが、K1にとってはcriticalではない

④余談:方法1

ちなみに、再び①の終了時の状態に戻して、K1の固定を外してバンドル調整を実施し、

K1を含めてバンドル調整

「歪曲プロット」のダイアログで、K1の推定値に加えて「誤差」(おそらく標準偏差の推定値)が確認できる。この値0.00779871が、真値(今回は0だが)や推定値-5.1518e-05に対してかなり大きいことから、K1は不定ではないにしろ、やはりよく言われている通り、平行光軸撮影ではK1の 推定のバラつきが大きいこと確認できる。なお「誤差」欄の値が本当に、分散共分散行列などから求めた推定値の標準偏差ならば、その確認は前記事の検証方法1「最適解における分散共分散行列において、θiの分散が非常に大きいかを調べる。」に相当する。

K1の推定値の標準偏差(おそらく)を確認

 

バンドル調整でパラメータの不定性を確認するにはどうするか

カメラパラメータをθ、それらの非線形関数である再投影誤差のRMSをg(θ)で表すと、バンドル調整は、g(θ)を最小化する非線形最小二乗法である。

本研究室ではバンドル調整の数値シミュレーションを多用するが、Critical Configurationによってθの要素の一部、例えばθiが不定となることがある。不定と疑われるθiが本当に不定かどうか、数値的に確認したいとき、どうすればよいのだろうか?

θの真値や最適解からθiだけを1次元に動かして、g(θ) VS θiの散布図を描くのは間違いだ。この方法では、「g(θ)が実はθiの関数ではなかった」というおっちょこちょいな場合を除き、不定性を見つけられないだろう。

なぜかと言えば、あるパラメータの不定性は、他のパラメータとの連動によって成り立つからだ。例えば平行光軸撮影によるfの不定性は、fの推定を誤っても、外部パラメータを連動して誤らせることでg(θ)を増やさずに済むから成り立っている。fだけ変えてg(θ) VS fの散布図を描けば、見事な谷が見えるだろう。

同じ意味で、g(θ)のθiによる二階偏微分が0かどうかを調べるのも間違いだ。θiの不定性とは言っても、θiを含む多数のパラメータが連動することで生じている不定性を、代表してθiの不定性とみなそうとしているわけなのだから、θiだけ見ていても真実は見えない。

 

ではどうすればよいのか?私が使えている方法は、次のようなものだ。

  1. 最適解における分散共分散行列において、θiの分散が非常に大きいかを調べる。
  2. 最適化を実施したときの目的関数値が、θiを誤った値に固定した場合とθiを正しい値に固定した場合で、有意に変わらないかを調べる。
  3. θiに誤った初期値を与えて最適化を実施したとき、明らかに誤った値に収束するかを調べる。

1はストレートだが、分散共分散行列を数値的に適切に求められるよう、方法に注意が必要だ。

2では、その他の誤差要因の影響がある限り、誤った値に固定した場合を何ケースが実施して、正しい場合より有意に大きいかどうかを判断する必要がある。

3で自信が持てるのは、その他の誤差要因や多峰性の影響を取り除ける場合に限られる。

 

2視点三角測量の精度の不思議

光軸が平行かつ、光軸と基線が直交するように撮影した2枚の写真を用いて、写真間で対応付けられた点Pの座標を推定する(2視線の交点の座標を計算する)「2視点三角測量」(two-view triangulation) の状況を考える。最も古典的な写真測量だ。

各写真における点Pの写った位置の観測誤差が、互いに独立で同じ分布(平均:0 かつ 標準偏差:σ )に従うと仮定したとき、光軸に平行な奥行き方向の分散は基線長に依存する。一方で、光軸と直交し撮像面に平行な2つの軸方向(撮像面の横方向、縦方向)の分散は基線長に依存しない

これらは、点Pの座標を求める式(2視線の交点の座標を表す式)に誤差伝搬の法則を適用することで導かれる。そのように理論的に導かれた分散は、Abdel-Aziz (1982)などに示されており、その平方根である標準偏差は次式の通りとなる。

  • 奥行き方向の座標推定の標準偏差=(奥行き方向の距離/基線長)× (奥行き方向の距離/焦点距離)×σ×√2
  • 撮像面の横方向・縦方向の座標推定の標準偏差= (奥行き方向の距離/焦点距離)×σ/√2

ただし、後者については本当は、点Pの撮像面の横方向・縦方向の座標に依存する。ここでは簡単のため、点Pの撮像面の横方向・縦方向の座標が基線の中点のそれらと一致する場合を考えている

 

基線長が短くなると、奥行き方向の精度が劣化することは、両目で遠くの物までの距離を掴みづらいことからも感覚的に理解できるし、「基線が短いと、両カメラの視線を作図したとき、その交点付近が『奥行き方向に長いX』になって、奥行き方向の位置の判読が難しくなるでしょう?」というイメージの説明でも一応納得がいく(本当は、観測誤差があるのは撮像面上なので、一様な太さのペンで描いた視線の太さを誤差とみなすのは厳密でないと思われる;最下部の図のように説明する方がよいはず)

 

しかし、基線長が短くなっても、奥行き方向でない方向の精度は変わらないというのは、実は私には不思議であることに最近気づいた。確かに理論的にはそう導けるのだが、感覚的にはしっくりきていなかった。数日考えてようやく作ったイメージの説明は、「基線が短いと、両カメラの視線を作図したとき、その交点が『奥行き方向に長いX』になって、奥行き方向の位置の判読に誤差が生じても、両視線の向きの違いが少ないわけだから、撮像面に平行な方向の座標の判断には影響がないでしょう?」というものだ。

 

なお、図を使えば、イメージをもっと明解に説明できる。下図では、真の投影位置pの左右に、投影位置を観測誤差εだけずらした2点"p-ε", "p+ε"を描き、それら3点を通る両カメラの視線の交点を使って、奥行き方向・撮像面方法の推定位置の変動幅を説明している。

2視点三角測量の精度の図解(基線が比較的長い場合)

2視点三角測量の精度の図解(基線が比較的長い場合)

 

2視点三角測量の精度の図解(基線が比較的短い場合)

2視点三角測量の精度の図解(基線が比較的短い場合)

今回のように基本的な内容でも、理解したつもりになってから何年も経って、理解が浅いことに気づいたりする。数式で導くだけで満足して、感覚的な理解を怠っていたことに気づくこともある。写真測量は3次元空間の話なのだから、本来は理論の代数的理解と、その幾何学的な意味の感覚的理解が、いつも車の両輪として伴っているべきなのだろう。

 

【研究メモ230427】Metashapeのバンドル調整での「追加の補正を調整」(Fit additional corrections)オプションとは

 

Metashapeでバンドル調整を行うためのダイアログ「カメラアラインメントを最適化」(旧「カメラを最適化」;英語ではOptimize Cameras)には、バージョン1.6.0で、「追加の補正を調整」(Fit additional corrections)というチェックボックスが実装されました。

 

Metashape 1.8.4の「カメラアラインメントを最適化」ダイアログ

その趣旨についてはこちらの記事で触れていましたが、ブラックボックスなのでその後2年以上、研究には用いていませんでした。

おそらく、高次の接線方向歪みのパラメータp3, p4をまとめたもの+αなのかなと想像していたところですが、どうも違うようだということが、今回判りました。

このオプションをONにした状態でバンドル調整を実施し、「カメラキャリブレーション」ダイアログから内部パラメータをxmlファイルにエクスポートすると、xmlファイルに

<corrections type="fourier">

というセクションが追加され、そこにたくさんの数字が並びます。fourierという文字列、そしてこのセクションの中身から、p3やp4ではなくフーリエ係数に見えます。

私の予想では、f, cx, k1, p1などの既存の他の内部パラメータをあてはめた(再投影誤差のRMSが最小になるように最適化した)ときの残差、つまり他の内部パラメータで説明しきれない再投影誤差)を、画像全体を1周期とみなしてフーリエ級数展開したときの、フーリエ係数を格納しているのではないかと思います。

内部パラメータをエクスポートしたXMLファイルの例

他の内部パラメータを係数とする「Brownのカメラモデル」は物理的な(光学的な)モデルであり、「画像保存時の前処理で歪み補正されてしまった画像」には上手くあてはまらない懸念がありますが (James et al., 2020)、

もし上記の解釈通り、「追加の補正」に使われるモデルがフーリエ級数展開ならば物理的なモデルではありませんので、「画像保存時の前処理で歪み補正されてしまった画像」にも特に問題なく適用できると考えられます。そのような画像を使うMetashapeプロジェクトでは、k1, k2, k3やp1, p2などは敢えて推定せず(アラインメント前から0に固定し、その後「カメラアラインメントを最適化」ダイアログでもチェックを外したままにしておき)、「追加の補正」だけにする方がよい場合もあり得るかと思います。

 

James, M. R., Antoniazza, G., Robson, S., and Lane, S. N. (2020) Mitigating systematic error in topographic models for geomorphic change detection: accuracy, precision and considerations beyond off-nadir imagery. Earth Surf. Process. Landforms, 45: 2251– 2271. https://doi.org/10.1002/esp.4878.