はじめに
天文学のデータ解析では、観測される天体の構造は PSF（Point Spread Function）によってぼやけます。したがって、観測データの広がりが PSF の影響によるものなのか、それとも天体の本来の構造によるものなのかを見極めることは、解析において重要な課題となります。
例えば、X線天文学では、PSF をシミュレーションや実測データをもとに見積もり、観測データの真の空間的スケールを推定します。その際、イメージデコンボリューション法を用いて直接真の天体像を推定することもありますが、ノイズの影響を受けやすいという課題があります。そのため、PSF の影響を直感的に理解し、頭の中で逆畳み込みをイメージできる ことが、より適切な解析手法の選択につながる場合があります。
本記事では、PSFと真の構造をガウシアン分布とみなし、観測される広がりの変化を数学的に整理し、可視化を通じて直感的に理解することを目指します。
本記事の実装コードは Google Colab からも確認できます。
観測される広がりの数式
天文学や画像処理において、観測されるデータは真の空間構造とPSFの畳み込みとして得られます。まず、畳み込みの数式を導入し、その後、PSFと真の構造をガウシアン分布と仮定した場合の広がりを考えます。
畳み込みとは
観測される画像 $I(x)$ は、真の構造 $f(x)$ とPSF $p(x)$ の畳み込みで表されます。
I(x) = (f * p)(x) = \int_{-\infty}^{\infty} f(t) p(x - t) dt
この畳み込みにより、PSFの影響で真の構造がぼやけます。
PSFと真の構造がガウシアンの場合
PSFと真の構造をガウシアン分布と仮定すると、それぞれ次のように表されます。
-
真の構造
f(x) = \frac{1}{\sqrt{2\pi} \sigma_{\text{true}}} \exp \left(-\frac{x^2}{2\sigma_{\text{true}}^2} \right) -
PSF
p(x) = \frac{1}{\sqrt{2\pi} \sigma_{\text{psf}}} \exp \left(-\frac{x^2}{2\sigma_{\text{psf}}^2} \right)
ガウシアン同士の畳み込みもガウシアンになり、観測される広がりの標準偏差は
\sigma_{\text{obs}} = \sqrt{\sigma_{\text{true}}^2 + \sigma_{\text{psf}}^2}
で表されます。
数式の可視化
$\sigma_{\text{obs}} = \sqrt{\sigma_{\text{true}}^2 + \sigma_{\text{psf}}^2}$の関係を可視化してみたいと思います。横軸に真の構造の標準偏差、縦軸に観測の標準偏差をとります。
数式可視化用のPythonコード
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 12
# パラメータ設定
sigma_true_array = np.linspace(0.01, 5, 100) # x は真の構造のスケール
sigma_psf = 1 # PSFの空間分解能
sigma_obs_array = np.sqrt(sigma_true_array**2 + sigma_psf**2)
# プロット
plt.figure(figsize=(7, 7))
plt.plot(sigma_true_array, sigma_obs_array, label=fr"With PSF: $\sigma_{{\text{{obs}}}} = \sqrt{{\sigma_{{\text{{true}}}}^2 + \sigma_{{\text{{psf}}}}^2}}$, ($\sigma_{{\text{{psf}}}}={sigma_psf}$)", color='k')
plt.plot(sigma_true_array, sigma_true_array, label=fr"Without PSF: $\sigma_{{\text{{obs}}}}=\sigma_{{\text{{true}}}}$", color='red', linestyle="--")
plt.title("Observed Standard Deviation with PSF")
plt.xlabel(r"$\sigma_{\text{true}}$")
plt.ylabel(r"$\sigma_{\text{obs}}$")
plt.legend()
plt.grid(True)
plt.gca().set_aspect('equal')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 12
# パラメータ設定
sigma_true_array = np.linspace(0.01, 5, 100) # x は真の構造のスケール
sigma_psf = 1 # PSFの空間分解能
sigma_obs_array = np.sqrt(sigma_true_array**2 + sigma_psf**2)
# 中心強度の計算
I_true_values = 1 / (np.sqrt(2 * np.pi) * sigma_true_array)
I_obs_values = 1 / (np.sqrt(2 * np.pi) * sigma_obs_array)
# プロット
plt.figure(figsize=(7, 7))
plt.plot(sigma_true_array, I_obs_values, label=fr"With PSF: $\sigma_{{\text{{obs}}}} = \sqrt{{\sigma_{{\text{{true}}}}^2 + \sigma_{{\text{{psf}}}}^2}}$, ($\sigma_{{\text{{psf}}}}={sigma_psf}$)", color='k')
plt.plot(sigma_true_array, I_true_values, label=fr"Without PSF: $\sigma_{{\text{{obs}}}}=\sigma_{{\text{{true}}}}$", color='red', linestyle="--")
plt.title("Observed Peak Intensity with PSF")
plt.xlabel(r"$\sigma_{\text{true}}$")
plt.ylabel(r"$I_{\text{obs}}(0)$")
plt.legend()
plt.grid(True)
plt.yscale('log')
plt.show()
出力結果
- PSFがある場合（下図は、$\sigma_{\text{psf}}=1$に固定して表示したもの）
- PSFがない場合
| 観測される広がりの変化 | 中心部の強度 |
|---|---|
 |
 |
この可視化から、以下のことが分かります。
- 真の構造がPSFより大きいと、PSFの影響は小さい
- 真の構造がPSFより小さいと、観測される広がりが大きくなる
実際のガウス分布で確認
ここでは、実際のガウス分布を用いて、PSFの影響を直感的に捉える練習をします。
PSFと真の空間構造のガウス分布の確認コード
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 12
# 1次元ガウシアン関数
def gaussian_1d(x, sigma):
return np.exp(-x**2 / (2 * sigma**2)) / (np.sqrt(2 * np.pi) * sigma)
# パラメータ設定
sigma_true = 3 # 真の構造の広がり
sigma_psf = 1 # PSFの広がり
sigma_obs = np.sqrt(sigma_true**2 + sigma_psf**2) # 観測される広がり
# x軸の範囲
x_values = np.linspace(-20, 20, 500)
# 各ガウシアンを計算
y_true = gaussian_1d(x_values, sigma_true)
y_true /= np.sum(y_true)
y_psf = gaussian_1d(x_values, sigma_psf)
y_psf /= np.sum(y_psf)
# 畳み込み演算
y_obs = np.convolve(y_true, y_psf, mode='same')
# プロット
plt.figure(figsize=(1.618*5, 5))
plt.plot(x_values, y_true, label=fr"$\sigma_{{\text{{true}}}} = {sigma_true}$", color='k')
plt.plot(x_values, y_psf, label=fr"$\sigma_{{\text{{psf}}}} = {sigma_psf}$", color='gray', linestyle=":")
plt.plot(x_values, y_obs, label=fr"$\sigma_{{\text{{obs}}}} = \sqrt{{\sigma_{{\text{{true}}}}^2 + \sigma_{{\text{{psf}}}}^2}}={sigma_obs:.2f}$", color='red', linestyle="--")
plt.xlabel("x")
plt.ylabel("Intensity")
plt.title("1D Gaussian Convolution: Effect of PSF")
plt.legend()
plt.grid(True)
plt.show()
- 真の空間スケール << PSFの空間分解能: $\sigma_{\text{psf}}=1, \sigma_{\text{true}}=0.5$
- 真の空間スケール = PSFの空間分解能: $\sigma_{\text{psf}}=1, \sigma_{\text{true}}=1$
- 真の空間スケール >> PSFの空間分解能: $\sigma_{\text{psf}}=1, \sigma_{\text{true}}=2$
パラメータを変えることで、PSFの影響を直感的に理解できます。実際に試してみると、ルートの項が効いてくるため、PSFと真の構造の広がりの関係が、直感的な予想と少し異なると感じる読者もいるかもしれません。
GIFアニメーションで確認
PSFの影響がどのように変化するかをアニメーションでも確認します。
GIFアニメーション作成コード
import imageio
import os
sigma_psf = 1
# sigma_true の範囲
sigma_true_values = np.arange(0.2, 5.2, 0.2) # 5 から 20 まで変化
# x軸の範囲
x_values = np.linspace(-20, 20, 500)
# 画像ファイルリスト
image_files = []
for sigma_true in sigma_true_values:
# 観測される広がりを計算
sigma_obs = np.sqrt(sigma_true**2 + sigma_psf**2)
# 各ガウシアンを計算
y_true = gaussian_1d(x_values, sigma_true)
y_true /= np.sum(y_true)
y_psf = gaussian_1d(x_values, sigma_psf)
y_psf /= np.sum(y_psf)
# 畳み込み演算
y_obs = np.convolve(y_true, y_psf, mode='same')
# プロット
plt.figure(figsize=(1.618 * 5, 5))
plt.plot(x_values, y_true, label=fr"$\sigma_{{\mathrm{{true}}}} = {sigma_true:.1f}$", color='k')
plt.plot(x_values, y_psf, label=fr"$\sigma_{{\mathrm{{psf}}}} = {sigma_psf}$", color='gray', linestyle=":")
plt.plot(x_values, y_obs, label=fr"$\sigma_{{\mathrm{{obs}}}} = \sqrt{{\sigma_{{\mathrm{{true}}}}^2 + \sigma_{{\mathrm{{psf}}}}^2}}={sigma_obs:.2f}$", color='red', linestyle="--")
plt.xlabel("x")
plt.ylabel("Intensity")
plt.title("1D Gaussian Convolution: Effect of PSF")
plt.legend()
plt.grid(True)
# 画像保存
image_path = os.path.join("./", f"frame_{sigma_true:.2f}.png")
plt.savefig(image_path)
image_files.append(image_path)
plt.close()
# GIF アニメーションを作成
imageio.mimsave("./gaussian_convolution.gif", [imageio.imread(img) for img in image_files], fps=3)
このGIFから、PSFによる広がりの変化を視覚的に確認できます。小さい構造ほどぼやけやすく、大きな構造では影響が小さいことが分かります。
まとめ
本記事では、PSF（Point Spread Function）が観測データに与える影響について、ガウシアン分布を用いて説明し、可視化を通じて直感的に理解できるようにしました。
結果として、小さい構造ほどPSFの影響を強く受け、十分に大きな構造では影響が小さくなる ことが確認できました。PSFの影響を正しく理解することで、観測データの解釈や分析の精度を高めることができます。
PSFの影響を意識することで、より正確なデータ解析が可能になります。本記事が、その理解の一助となれば幸いです。