佐藤の積分公式

数値積分における佐藤の積分公式(Sato integrator)とは、重みつきの直交関数系フィルターがある場合、それを通常の長方形公式や台形公式で数値計算すると、バイアスが生じる場合があるので、そのようなバイアスが生じないように設計されたインテグレーター(anti-bias integrator)である。

例えば、1次元のガウシアン型の関数

$g(m,\rho ,x)=\frac{N}{\sqrt{2\pi \rho }}\exp (-\frac{(x-m{)}^2}{2{\rho }^2})$

の形をした分布関数が、1次元の格子上に落ちたとしよう。 格子のサイズを1として、その格子は、その区間上の積分

$I(x)=\underset{x}{\overset{x+1}{\int }}g(m,\rho ,t)dt$

の値を持っている。 実際の画像データの場合は、これにフォトンノイズやホワイトノイズが加わる。 フォトンノイズは、$I(x)$の平方根に比例する強度のノイズ、ホワイトノイズは一定の強度のノイズであるため、ウィンドウを大きくし過ぎるとノイズが強くなって測定精度を低下されることになり、適当な大きさのウィンドウをかける必要がある。

このデジタル化された分布関数に対して、ガウシアン型の重みを持った直交関数系フィルターをかけて、そのn次のモーメントを求める。 その時、そのn次関数フィルター(正規直交関数系)は、エルミート関数(Hermite functions)で次のように表される。

$s_n(x)=\sqrt{\frac{\sqrt{2+w^2}}{n!w}}{H}_n\left(\frac{\sqrt{2+w^2}}{w}\right)\exp (-\frac{x^2}{2w^2})$

ここで、$w$は、ガウシアン型フィルターの幅$\sigma$とイメージサイズ$\rho$の比で、$s_n(x)$は次の直交関係を満たす。

$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{s}_n(x)\overline{s_m(x)}g(0,1,x)dx={\delta }_{mn}$

任意の次数のフィルターは、$s_n(x)$の線形結合で作られる。

まず、0次のフィルターを考えると、それは測定する像分布の強度(フォトン数)を測定するフィルターである。 その推定分散を最小にする関数は、

$f_0(x)=\frac{1+w^2}{w}\exp (-\frac{x^2}{2w^2{\rho }^2})$

で与えられ、この時の推定分散は、次のように評価される。

${\sigma }_N^2=\frac{(1+w^2)N}{w\sqrt{2+w^2}}(1+\sqrt{\pi (2+w^2)}\frac{B\rho }{N})$

ここで、第1項はフォトンノイズによる項、第2項はバックグラウンドのホワイトノイズ$B$による項である。 また、推定値には、次のバイアス

${\mathrm{\Delta }}_N=-\frac{w^3\sqrt{2+w^2}}{12(1+w^2{)}^2}$

が生じる。 このように、0次の(一般的には偶数次の)フィルターでは測定値にバイアスが生ずる。 観測値はピクセルでデジタル化されているため、それに掛けるフィルターの値を調整することでしか、このバイアスを消去できない。 このバイアスが生じないようにするインテグレーターが佐藤公式である。

有限ステップのインテクレーターには、長方形公式、台形公式、シンプソン公式などがよく知られているが、佐藤公式の特徴は、隣り合う2区間(2ピクセル)をワンセットで考えることである。 その公式は、長方形公式、台形公式、シンプソン公式をそれぞれ

$\left[\begin{array}{ccc}0& 1& 0\end{array}\right]$

$\left[\begin{array}{ccc}1& 0& 1\end{array}\right]$

$\left[\begin{array}{ccc}1& 4& 1\end{array}\right]$

とすると、

$\left[\begin{array}{ccc}0& 4& -1\end{array}\right]\left[\begin{array}{ccc}-1& 4& 0\end{array}\right]$

である。 2次元では、

$\left[\begin{array}{ccc}0& 0& 0\\ 0& 16& -4\\ 0& -4& 1\\ \end{array}\right]\left[\begin{array}{ccc}0& 0& 0\\ -4& 16& 0\\ 1& -4& 0\\ \end{array}\right]$

$\left[\begin{array}{ccc}0& -4& 1\\ 0& 16& -4\\ 0& 0& 0\end{array}\right]\left[\begin{array}{ccc}1& -4& 0\\ -4& 16& 0\\ 0& 0& 0\end{array}\right]$

となる。 この時、n次のモーメントの推定値は、(1次元の場合)

$M_n={\mathrm{\Sigma }}_x(\frac{4s_n(x+1/2)-s_n(x+1)}{3}I(x)+\frac{-s_n(x+1)+4s_n(x+3/2)}{3}I(x+1))$

となる。 0次と1次のモーメントの比より、中心位置は、$x_0=M_1/M_0$と求められる。

この公式は、星像のように点光源に近い天体像のアストロメトリー(位置、明るさ、形状の測定)に向いているが、裾野が減衰する一般の分布関数にも精度よく適用できる。

出典

• Filtering theory with orthogonal functions for two-dimensional astrometry, Sato I. & Yoshizawa M. (1995), Astronomy & Astrophysical Journal Supplement Series, 113, pp.185-197

関連項目

• 長方形公式
• 台形公式
• シンプソンの公式
• 数値積分
• アストロメトリー