Statistics of Centroids
位置天文学ではイメージセンサで捉えた星像から星の中心位置を正確に求めることが必要になることがあります. ここでは中心位置 (重心) の推定精度について解析的にまとめてみます.
データと問題設定
ここでは観測で得られるデータを抽象化して考えます. 実際のデータは 2 次元画像ですが, 空間方向に相関はないと考えてまずは 1 次元で考えます. 得られる情報はある位置 \(x\) に入射した光子数 (電子数) \(N_x\) のペア \(\{x, N_x\}\) とします.
求めたい値は星像の中心位置 (重心) \(\bar{x}\) です. 以下の式で定義されます.
\[
\bar{x} = \frac{\sum_x x N_x}{\sum_x N_x}.
\]
ケース 0: 天体のみ
まずは \(N_x\) のソースとして天体からの光のみを考えます. 観測装置が天体から受けとる光子の総数の期待値を \(N_0\) と置きます. この光子が位置 \(x\) に分配される割合を \(\eta(x)\) と書きます. ただし, 光子数は Poisson 統計に従って揺らぎます. ゆらぎを表すための項として \(\varepsilon_x\) を加えます.
\[
N_x = N_0 \eta(x) + \varepsilon_x.
\]
\(\eta(x)\) は分配関数なので総和は 1 になるように規格化されます. また \(\eta(x)\) には星像の形状の情報が含まれています. ここでは 1 次と 2 次のモーメントが以下のようになると仮定します. ここで座標系の原点は星の位置からさほど離れていないとして \(\mu \ll \sigma\) という状況を想定しておきます.
\[
\begin{aligned}
\int\mathrm{d}x\, \eta(x) & \simeq \sum_x \eta(x) = 1 \\
\int\mathrm{d}x\, x\eta(x) & \simeq \sum_x \eta(x) = \mu \\
\int\mathrm{d}x\, x^2\eta(x) & \simeq \sum_x \eta(x) = \sigma^2 + \mu^2
\end{aligned}
\]
星のプロファイルの模式図
ゆらぎについては Poisson 統計にしたがうよう以下のように定義します. ここで \(\langle \cdot \rangle\) は繰り返し観測をおこなったときの \(\cdot\) の期待値を表しています.
\[
\begin{aligned}
\langle \varepsilon_x \rangle & = 0 \\
\langle \varepsilon_x^2 \rangle & = N_0\eta(x) \\
\langle \varepsilon_x\varepsilon_{x'} \rangle & = 0 \quad (x \neq x')
\end{aligned}
\]
以上の問題設定で中心位置 \(\bar{x}\) の期待値と分散を考えていきます.
期待値の評価
期待値は以下のように計算できます.
\[
\langle \bar{x} \rangle =
\left\langle
\frac{\sum_x x N_x}{\sum_x N_x}
\right\rangle
\]
ここでさっそくつまづきます. 確率的に揺らがないものを \(\langle \cdot \rangle\) の外に出したいのですが, 分母にも分子にも確率的に揺らぐ成分 \(\varepsilon_x\) が含まれているためにややこしいことになっています. 仮に分母・分子をそれぞれの期待値を撮ることができれば良いのですが, 残念ながら成立しません.
\[
\langle \bar{x} \rangle =
\left\langle
\frac{\sum_x x N_x}{\sum_x N_x}
\right\rangle
\neq
\frac{\left\langle \sum_x x N_x \right\rangle}
{\left\langle \sum_x N_x \right\rangle}
\]
成立しない理由は分母と分子の相関です. 分母を \(X\), 分子を \(Y\) と置いて整理します.
\[
\begin{aligned}
X & = \sum_x N_x
= \sum_x N_0 \eta(x) + \sum_x \varepsilon_x
= N_0 + \sum_x \varepsilon_x \\
Y & = \sum_x x N_x
= \sum_x N_0 x \eta(x) + \sum_x x \varepsilon_x
= \mu N_0 + \sum_x x \varepsilon_x
\end{aligned}
\]
定数は気にせずゆらぎの項だけに注目します. ゆらぎの項 \(\varepsilon_x\) は \(x \simeq \mu\) で大きな分散値を持ちます. \(X\), \(Y\) のゆらぎ成分の大半は \(x \simeq \mu\) の項が担っていると言えます. そこでゆらぎ成分だけに注目すると,
\[
\delta X \simeq \varepsilon_\mu, \quad
\delta Y \simeq \mu\varepsilon_\mu
\]
となるので \(X\), \(Y\) には明らかな相関が存在します. このため, 積の期待値をそれぞれの期待値の積で計算することができません.
さて, 天下り的ですが中心位置 \(\bar{x}\) ではなく \(\bar{x} - \mu\) の期待値を考えててみましょう.
\[
\langle \bar{x} - \mu \rangle =
\left\langle
\frac{\sum_x (x{-}\mu) N_x}{\sum_x N_x}
\right\rangle
\]
同様に分母と分子を \(X\), \(Y\) と置きます.
\[
X = N_0 + \sum_x \varepsilon_x, \quad
Y = \mu N_0 + \sum_x (x{-}\mu) \varepsilon_x
\]
\(X\), \(Y\) の共分散を計算します.
\[
\begin{aligned}
\mathrm{Cov}(X, Y)
& = \bigl\langle
(X - \langle X \rangle)
(Y - \langle Y \rangle)
\bigr\rangle \\
& = \left\langle
\sum_x \varepsilon_x
\sum_x (x{-}\mu) \varepsilon_x
\right\rangle \\
& = \left\langle
\sum_x (x{-}\mu) \varepsilon_x^2
+
2 \sum_{x \neq x'} (x{-}\mu)\varepsilon_x\varepsilon_{x'}
\right\rangle \\
& =
\sum_x (x{-}\mu) \bigl\langle \varepsilon_x^2 \bigr\rangle
+ 2 \sum_{x \neq x'} (x{-}\mu)
\bigl\langle\varepsilon_x\varepsilon_{x'} \bigr\rangle\\
& =
\sum_x (x{-}\mu) N\eta(x) \\
& = 0
\end{aligned}
\]
すると, 上記のようにゼロとなり \(X\), \(Y\) が無相関だとみなせることがわかります. \(X\) の変動成分はおもに \(x \simeq \mu\) 付近の項に由来しますが, \(Y\) の変動成分は \(x = \mu\) でゼロなので実質的に相関がなくなったように振る舞うことになります.
そこで, この性質を利用して以下のように近似します.
\[
\begin{aligned}
\langle \bar{x} - \mu \rangle
& = \left\langle
\frac{\sum_x (x{-}\mu) N_x}{\sum_x N_x}
\right\rangle
\simeq
\left\langle \frac{\sum_x (x{-}\mu) N_x}{1} \right\rangle
\left\langle \frac{1}{\sum_x N_x} \right\rangle
\simeq
\frac{\left\langle \sum_x (x{-}\mu) N_x \right\rangle}
{\left\langle \sum_x N_x \right\rangle}
\end{aligned}
\]
最後の近似は \(\sum_x N_x\) のゆらぎが平均値に対して十分に小さいときに, 逆数の期待値と期待値の逆数がほぼ一致することを使用しています. それぞれ計算すると以下のとおりです.
\[
\begin{aligned}
\left\langle \sum_x (x{-}\mu) N_x \right\rangle
& = N_0 \sum_x (x{-}\mu)\eta(x)
+ \sum_x (x{-}\mu){\langle \varepsilon_x \rangle}
= N_0 (\mu - \mu) = 0 \\
\left\langle \sum_x N_x \right\rangle
& = \sum_x N_0 \eta(x) + \sum_x {\langle \varepsilon_x \rangle}
= N_0
\end{aligned}
\]
これより \(\langle \bar{x} - \mu \rangle = 0\) となり, 以下の式を得ます.
\[
\langle \bar{x} \rangle
= \left\langle
\frac{\sum_x x N_x}{\sum_x N_x}
\right\rangle
= \mu
\]
定義した計算で星像の中心位置 \(\mu\) をバイアスなく推定できることがわかります.
分散の評価
次に \(\bar{x}\) の分散 \(\mathrm{Var}(\bar{x})\) を評価します. 先ほど \(\langle\bar{x}\rangle\) ではなく\(\langle \bar{x} - \mu \rangle\) を評価したほうが統計的な性質が扱いやすくなることを確認しました. 数を引いても分散の値は変化しないことから, ここでも \(\bar{x} - \mu\) の分散を評価することにします.
\[
\begin{aligned}
\mathrm{Var}\left(\bar{x} - \mu \right)
& = \mathrm{Var}\left(
\frac{\sum_x (x{-}\mu) N_x}{\sum_x N_x}
\right) \\
& = \left\langle
\frac{\left( \sum_x (x{-}\mu) N_x \right)^2}{\left( \sum_x N_x \right)^2}
\right\rangle
- \underset{\qquad =0}{\underline{\left\langle
\frac{ \sum_x (x{-}\mu) N_x }{ \sum_x N_x }
\right\rangle}}^2 \\
& = \left\langle
\left( \frac{\sum_x (x{-}\mu) N_x}{1} \right)^2
\right\rangle
\left\langle
\left( \frac{1}{\sum_x N_x} \right)^2
\right\rangle \\
& \simeq \frac{1}{N_0^2}
\left\langle
\left( \sum_x (x{-}\mu) N_x \right)^2
\right\rangle \\
& = \frac{1}{N_0^2}
\mathrm{Var}\left( \sum_x (\bar{x} - \mu)N_x \right) \\
& = \frac{1}{N_0^2}
\sum_x (x{-}\mu)^2 \mathrm{Var}(N_x)
\end{aligned}
\]
念のために数値実験で上記の式が成立することを確認しました. 数値実験の流れは以下のとおりです. 正規分布にしたがうようなプロファイル \(\eta(x)\) を作成して 10000 個の光子を分配します. Poisson 分布にしたがう乱数でフラックスに対応した光子をサンプリングする操作を 10000 回繰り返しました. このデータを使用して上記の操作の妥当性を確認してみます.
pkg load statistics
N0 = 1e4;
mu = 2.5
sig = 1.00
x = [-10:10];
profile = @(mu, sig) normpdf(x, mu, sig);
Nx=[];
for n=1:10000;
Nx(end+1,:) = poissrnd(N0 * profile(mu, sig));
endfor
\(\mathrm{Var}(\bar{x})\) をそのまま評価した結果は以下のとおりです. この値はおおよそ \(\sigma^2 / N_0\) に対応します.
octave> var(sum(Nx .* x) ./ sum(Nx, 2), 2)
ans = 9.9715e-05
\(\mathrm{Var}(\bar{x})\) の分子を \(N_0\) に置き換えて評価すると値が大きくなります. つまり勝手に \(N_0\) に置き換えてはいけません. なお, この値はおおよそ \((\mu^2 + \sigma^2)/N_0\) に対応します.
octave> var(sum(Nx .* x) ./ N0, 2)
ans = 7.2414e-04
しかし, \(\mathrm{Var}(\bar{x})\) の分子を \(N_0\) に置き換えて評価した場合には直接計算した結果とほぼ同じ値が返ってきます. 相関を排除したことで分母と分子でつじつま合わせをする必要がなくなったため, 別々に計算しても問題なくなったと捉えています.
octave> var(sum((Nx .* (x - mu)) ./ N0, 2))
ans = 9.9693e-05
さらに, この状態から分散の計算と和の計算を入れ替えてもほぼ同じ値が得られるとわかります. 2 番目の式は分散の定義に従って整理しただけですが, 分散を評価するときには \(N_x\) だけを考えればいいことがわかります.
octave> sum(var((Nx .* (x - mu)) ./ N0), 2)
ans = 9.9751e-05
octave> sum(var(Nx) .* (x - mu).^2, 2) / N0**2
ans = 9.9751e-05
以上を整理すると \(\bar{x} - \mu\) の分散について以下の近似が成り立ちます.
\[
\mathrm{Var}\left(\bar{x} - \mu \right)
\simeq \frac{1}{N_0^2}
\sum_x (x{-}\mu)^2 \mathrm{Var}\left(N_x\right)
\]
\(\mathrm{Var}(N_x) = \mathrm{Var}(\varepsilon_x) = N_0\eta(x)\) なので, 以下の式を得ます.
\[
\mathrm{Var}(\bar{x})
= \mathrm{Var}\left(\bar{x} - \mu \right)
\simeq \frac{1}{N_0} \sum_x (x{-}\mu)^2 \eta(x)
= \frac{\sigma^2}{N_0}.
\]
光子数 \(N_0\) の数に応じて中心位置推定精度が向上することが示せました.
ケース 1: 無相関なノイズ
先ほどは光子が Poisson 統計にしたがうことに依る影響のみを考えていました. ここでは天体からの光子とは関係なく場所によってほぼ一様に生じるノイズを考えます. これは検出器の読み出しノイズ, ダークカレントのゆらぎ, 背景光 (迷光) のゆらぎなどに対応します.
\[
N_x = N_0 \eta(x) + \varepsilon_x + \xi_x
\]
このノイズは以下の性質を持つとします.
\[
\begin{aligned}
\langle \xi_x \rangle & = 0 \\
\langle \xi_x^2 \rangle & = \delta^2 \\
\langle \xi_x\xi_{x'} \rangle & = 0 \quad (x \neq x')
\end{aligned}
\]
ノイズの分散は場所に依らず一定値 \(\delta^2\) をとります. また, ノイズ同士の相関も存在しない (\(\langle \varepsilon_x \xi_x \rangle = 0\)) と仮定します. このノイズの影響下で中心位置の推定精度がどのように変わるかを定式化します.
期待値の評価
ケース 0 と同様の近似をもちいて \(\bar{x} - \mu\) の平均値を評価します.
\[
\begin{aligned}
\langle \bar{x} - \mu \rangle
& \simeq
\frac{\left\langle \sum_x (x{-}\mu) N_x \right\rangle}
{\left\langle \sum_x N_x \right\rangle} \\
& = \frac{1}{N_0}
\left\langle
\sum_x (x{-}\mu)\left(
N_0 \eta(x) + \varepsilon_x + \xi_x
\right)
\right\rangle \\
& = \frac{1}{N_0}
\left(
\sum_x (x{-}\mu)
\left( N_0 \eta(x)
+ \underset{=0}{\underline{\langle\varepsilon_x\rangle}}
+ \underset{=0}{\underline{\langle\xi_x\rangle}}
\right)
\right) = 0
\end{aligned}
\]
ノイズの影響があっても期待値は変わらないことが確認できました.
分散の評価
分散もケース 0 と同様の近似をもちいて評価します.
\[
\begin{aligned}
\mathrm{Var}\left(\bar{x} - \mu \right)
& \simeq \frac{1}{N_0^2}
\sum_x (x{-}\mu)^2 \mathrm{Var}\left(N_x\right) \\
& = \frac{1}{N_0^2}
\sum_x (x{-}\mu)^2 \mathrm{Var}\bigl(
N_0 \eta(x) + \varepsilon_x + \xi_x
\bigr) \\
& = \frac{1}{N_0^2}
\sum_x (x{-}\mu)^2 \bigl(
\mathrm{Var}(\underset{\text{= const.}}{\underline{N_0 \eta(x)}})
+ \mathrm{Var}(\varepsilon_x)
+ \mathrm{Var}(\xi_x)
\bigr) \\
& = \frac{1}{N_0^2}
\sum_x (x{-}\mu)^2 \bigl(
N_0\eta(x) + \delta^2
\bigr) \\
& = \frac{\sigma^2}{N_0}
+ \frac{\delta^2}{N_0^2}\sum_x (x{-}\mu)^2
\end{aligned}
\]
最後の項 \(\sum_x(x{-}\mu)^2\) の大きさは足し合わせる \(x\) の範囲に依存します. ここでは, 解析対象が星像の大きさで決まると仮定します. \(x \in [\mu{-}\alpha\sigma, \mu{+}\alpha\sigma]\) として, 和を積分に置き換えることで以下のように見積もりました.
\[
\sum_x(x{-}\mu)^2
\simeq \int_{\mu{-}\alpha\sigma}^{\mu{+}\alpha\sigma}
\!\!\!\!\!\!\mathrm{d}x \, (x-\mu)^2
= \frac{2}{3}(\alpha\sigma)^3
\]
以上をまとめて星像中心位置の分散は以下のようになります.
\[
\mathrm{Var}\left(\bar{x}\right)
\simeq \frac{\sigma^2}{N_0}
+ \frac{2}{3}\frac{\delta^2}{N_0^2}(\alpha\sigma)^3
\]
新しく増えた項がノイズに依る影響です. 星からの光子数 \(N_0\) が多くなると急速に影響が小さくなりますが, 星像の大きさ (解析に使用する領域) が大きくなるとノイズの影響が強く出ます. 画像を広くとればとるほど多くのノイズをかき集めてしまうことに加えて, 重心計算ではノイズに距離を掛け算した量が考慮されるために \(\sigma^3\) の依存性を持つと解釈できます.
ケース 2: フラット補正誤差
イメージセンサの各ピクセルはやってきた光子を 100 % の効率で電子に変換することが理想ですが, 実際にはピクセルごとに個性があります. この個性はフラットフレーム補正をすることによって修正することができますが, 完全に補正しきれるわけではありません. ここではフラット補正処理の残存誤差によって効率を間違えて解析したケースを想定します.
\[
N_x = N_0 \eta(x) + \varepsilon_x + \phi_x
\]
新たに加えられたノイズ \(\phi_x\) の性質を考えます. フラット補正の残存誤差によって効率を \(1 + f_x\) 倍間違えたと表現します. 一般に \(f_x \ll 1\) です. \(f_x\) は位置 \(x\) に固定されているので同じ観測をしたときに常に同じ値が返ってきます. これまでのケースと同じように定義すると以下のような性質を持ちます.
\[
\langle \phi_x \rangle \simeq f_x N_0\eta(x), \quad
\mathrm{Var}(\phi_x) \simeq f_x^2 N_0\eta(x)
\]
\(\langle \phi_x \rangle \neq 0\) なので期待値の推定に固定のバイアスを与えます. 一方で \(f_x \ll 1\) なので分散にはほとんど寄与しません. フラット補正誤差が与える影響としては正しいのですが, 位置推定の精度に与える影響を見るという点では少し扱いづらくなります.
実際の観測では「測定を繰り返す」という作業はかならずしも検出器上の同一のピクセルで観測することを意味しません. そこで, \(\phi_x\) に対してのみ \(\langle \cdot \rangle\) の定義を「繰り返し観測をおこなったときの \(\cdot\) の期待値 (ただし毎回イメージセンサの異なるピクセルを使用したとする)」と変更します.
あらためて \(\phi_x\) の性質を以下のように定義します.
\[
\begin{aligned}
\langle \phi_x \rangle & = 0 \\
\langle \phi_x^2 \rangle & = (\kappa N_0 \eta(x))^2 \\
\langle \phi_x\phi_{x'} \rangle & = 0 \quad (x \neq x')
\end{aligned}
\]
\(\kappa\) はフラット補正誤差の典型的な大きさを表す量です. \(\kappa\) をもちいて \(f_x\) を定義すると \(f_x \sim \mathcal{N}(0, \kappa)\) となります. これまでと同様にノイズ間の相関は考えません.
期待値の評価
ケース 0 と同様の近似をもちいて \(\bar{x} - \mu\) の平均値を評価します. これまで同様にノイズの影響があっても期待値は変わりません.
\[
\begin{aligned}
\langle \bar{x} - \mu \rangle
& \simeq
\frac{\left\langle \sum_x (x{-}\mu) N_x \right\rangle}
{\left\langle \sum_x N_x \right\rangle} \\
& = \frac{1}{N_0}
\left\langle
\sum_x (x{-}\mu)\left(
N_0 \eta(x) + \varepsilon_x + \phi_x
\right)
\right\rangle \\
& = \frac{1}{N_0}
\left(
\sum_x (x{-}\mu)
\left( N_0 \eta(x)
+ \underset{=0}{\underline{\langle\varepsilon_x\rangle}}
+ \underset{=0}{\underline{\langle\phi_x\rangle}}
\right)
\right) = 0
\end{aligned}
\]
分散の評価
分散も同様に評価します.
\[
\begin{aligned}
\mathrm{Var}\left(\bar{x} - \mu \right)
& \simeq \frac{1}{N_0^2}
\sum_x (x{-}\mu)^2 \mathrm{Var}\left(N_x\right) \\
& = \frac{1}{N_0^2}
\sum_x (x{-}\mu)^2 \mathrm{Var}\bigl(
N_0 \eta(x) + \varepsilon_x + \phi_x
\bigr) \\
& = \frac{1}{N_0^2}
\sum_x (x{-}\mu)^2 \bigl(
\mathrm{Var}(\underset{\text{= const.}}{\underline{N_0 \eta(x)}})
+ \mathrm{Var}(\varepsilon_x)
+ \mathrm{Var}(\phi_x)
\bigr) \\
& = \frac{1}{N_0^2}
\sum_x (x{-}\mu)^2 \left(
N_0\eta(x) + (\kappa N_0 \eta(x))^2
\right) \\
& = \frac{\sigma^2}{N_0}
+ {\kappa^2}\sum_x (x{-}\mu)^2\eta(x)^2
\end{aligned}
\]
最後の項をおおざっぱに評価するため \(\eta(x)\) の関数系を中心が \(\mu\) で標準偏差が \(\sigma\) の Gauss 分布 \(\mathcal{N}(x; \mu, \sigma)\) で表現されると仮定して \(\eta(x)^2\) を計算します.
\[
\begin{aligned}
\eta(x)^2
& = \left[ \frac{1}{\sigma\sqrt{2\pi}}
\exp\left(
-\frac{1}{2}\left( \frac{x-\mu}{\sigma} \right)^2
\right) \right]^2 \\
& = \frac{1}{2\sigma\sqrt{\pi}} \frac{1}{(\sigma/\sqrt{2})\sqrt{2\pi}}
\exp\left(
-\frac{1}{2}\left( \frac{x-\mu}{\sigma/\sqrt{2}} \right)^2
\right) \\
& = \frac{1}{2\sigma\sqrt{\pi}}
\mathcal{N}\left(x; \mu, \frac{\sigma}{\sqrt{2}}\right)
\end{aligned}
\]
これより最終項は以下のようになります.
\[
{\kappa^2}\sum_x (x{-}\mu)^2\eta(x)^2
= \frac{\kappa^2}{2\sigma\sqrt{\pi}}
\sum_x (x-\mu)^2 \mathcal{N}\left(x; \mu, \frac{\sigma}{\sqrt{2}}\right)
= \frac{\kappa^2\sigma}{4\sqrt{\pi}}
\]
以上をまとめて星像中心位置の分散は以下のようになります.
\[
\mathrm{Var}\left(\bar{x}\right)
\simeq \frac{\sigma^2}{N_0}
+ \frac{\kappa^2\sigma}{4\sqrt{\pi}}
\]
フラット補正の残存誤差に依る影響は光子の総数 \(N_0\) には依存しません. これは光子数が多いほどノイズの影響が比例して増えるため, ノイズに対する S/N 比が変わらないためだと解釈できます. また, 星像が大きくなるとフラット残存誤差による影響は大きくなります. 中心位置から離れた位置にあるピクセルの影響を拾いやすくなるためだと考えられます.
まとめ
以上の結果をまとめます. 天体由来の Poisson ノイズ, 観測領域に一様に発生するノイズ, フラット補正の残存誤差を考慮に入れたケースを考えます. 位置 \(x\) で検出した光子数 (電子数) を以下のように表現します. それぞれのノイズの性質はこれまで定義したとおりです.
\[
N_x = N_0 \eta(x) + \varepsilon_x + \xi_x + \phi_x
\]
このデータから重心 \(\bar{x}\) を計算することで天体の中心位置を推定します. 推定した位置の信頼性をを推定値の分散 \(\mathrm{Var}(\bar{x})\) で評価します. ノイズ間の相関がないと思うと各ノイズ源が中心位置推定精度に与える影響は以下の式で評価できます.
\[
\mathrm{Var}(\bar{x}) \simeq
\frac{\kappa^2\sigma}{4\sqrt{\pi}}
+ \frac{\sigma^2}{N_0}
+ \frac{2}{3}\frac{\delta^2(\alpha\sigma)^3}{N_0^2}
\]
この式はそれぞれのノイズ源が推定精度にどのような影響を与えるのかを定性的に評価するためのものです. 実際の測定との定量的な評価には使えませんが, 測定精度がどのような関数系にしたがうかの根拠にはなると思います. 関数系だけに注目して整理すると以下のようになります.
\[
\mathrm{Var}(\bar{x}) \simeq
A_0 \kappa^2\sigma
+ A_1 \frac{\sigma^2}{N_0}
+ A_2 \frac{\delta^2\sigma^3}{N_0^2}
\]
2 次元イメージセンサへの拡張
最後に 2 次元イメージセンサで評価したときの影響を考えてみます. これまでの式を 2 次元イメージセンサに拡張するために, 光子数 (電子数) の分布を 2 次元で考えます.
\[
N_{x,y} = N_0 \eta(x,y) + \varepsilon_{x,y} + \xi_{x,y} + \phi_{x,y}
\]
2 次元にしてもそれぞれのノイズが従う関係は変わりません.
\[
\begin{aligned}
\langle \varepsilon_{x,y}^2 \rangle & = N_0\eta(x,y) \\
\langle \xi_{x,y}^2 \rangle & = \delta^2 \\
\langle \phi_{x,y}^2 \rangle & = (\kappa N_0 \eta(x,y))^2
\end{aligned}
\]
星像の中心位置 (重心) \(\bar{x}\) には \(y\) 方向の足し合わせも加わります.
\[
\bar{x} = \frac{\sum_x\sum_y x N_{x,y}}{\sum_x\sum_y N_{x,y}}.
\]
2 次元の分配関数 \(\eta(x,y)\) は \(\eta(x)\) と \(\eta(y)\) の積に分解できると仮定します. 以上の関係をもちいて 2 次元イメージセンサの場合の中心位置推定精度を考えます. \(\varepsilon_{x,y}\) に関係する項は以下のように整理でき, まったく同じ結果になります.
\[
\begin{aligned}
& \sum_x (x{-}\mu)^2 \sum_y \langle \varepsilon_{x,y}^2 \rangle \\
& \qquad\qquad = \sum_x (x{-}\mu)^2 \sum_y N_0 \eta(x,y) \\
& \qquad\qquad = N_0 \sum_x (x{-}\mu)^2 \eta(x) \\
& \qquad\qquad = N_0 \sigma^2
\end{aligned}
\]
次に \(\xi_{x,y}\) に関係する項を計算します. \(\sum_x\sum_y(x{-}\mu)^2\) の大きさは足し合わせる \(x\), \(y\) の範囲に依存します. ここでは 1 次元での検討と同様に解析対象が星像の大きさで決まると仮定して \(x \in [\mu{-}\alpha\sigma, \mu{+}\alpha\sigma]\), \(y \in [{-}\alpha\sigma, {+}\alpha\sigma]\) の範囲での積分から見積もりました.
\[
\begin{aligned}
& \sum_x (x{-}\mu)^2 \sum_y \langle \xi_{x,y}^2 \rangle \\
& \qquad\qquad = \delta^2 \sum_x \sum_y (x{-}\mu)^2 \\
& \qquad\qquad \simeq \delta^2
\int^{\mu{+}\alpha\sigma}_{\mu{-}\alpha\sigma}
\!\!\!\!\!\mathrm{d}x (x{-}\mu)^2
\int^{{+}\alpha\sigma}_{{-}\alpha\sigma}
\!\!\!\!\!\mathrm{d}y \\
& \qquad\qquad \simeq \frac{4}{3}(\alpha\sigma)^4
\end{aligned}
\]
このノイズはイメージセンサのどこでも同じように影響します. 解析する領域が現在考えている軸と直行する方向に広がると, その分だけ多くのノイズを拾うことになります. 1 次元のときと比べて \(2\alpha\sigma\) だけ変わっているのは面積で計算していることと関係があります.
最後に \(\phi_{x,y}\) に関係する項を計算します.
\[
\begin{aligned}
& \sum_x (x-\mu)^2 \sum_y \langle \phi_{x,y}^2 \rangle \\
& \qquad\qquad = \sum_x (x-\mu)^2 \sum_y \kappa^2N_0^2\eta(x,y)^2 \\
& \qquad\qquad \simeq \sum_x (x-\mu)^2 \sum_y
\frac{\kappa^2N_0^2}{4\pi\sigma^2} \eta'(x)\eta'(y) \\
& \qquad\qquad = \frac{\kappa^2N_0^2}{4\pi\sigma^2} \sum_x (x-\mu)^2 \eta'(x) \\
& \qquad\qquad = \frac{\kappa^2N_0^2}{8\pi}
\end{aligned}
\]
フラット残存誤差のときに使用した近似 (\(\eta(x) \sim \mathcal{N}(\mu, \sigma)\)) をここでも採用して計算しました. ここで \(\eta'(x)\) は標準偏差が \(\sigma/\sqrt{2}\) になった分配関数を意味しています.
\[
\sum_y \eta(x,y)^2
= \left(\sum_y\frac{1}{2\sigma\sqrt{\pi}}\right)^2 \eta'(x,y)
=
\]
これより 2 次元イメージセンサでの解析では以下の関係式が導かれます.
\[
\mathrm{Var}(\bar{x}) \simeq
\frac{\kappa^2}{8\sqrt{\pi}}
+ \frac{\sigma^2}{N_0}
+ \frac{4}{3}\frac{\delta^2(\alpha\sigma)^4}{N_0^2}
\]
さきほどと同様に各ノイズ源の依存性だけに注目して整理すると以下のようになります.
\[
\mathrm{Var}(\bar{x}) \simeq
A_0 \kappa^2
+ A_1 \frac{\sigma^2}{N_0}
+ A_2 \frac{\delta^2\sigma^4}{N_0^2}
\]
Last update:
November 25, 2023