はじめに

天文学のデータ解析では、観測される天体の構造は 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コード

PSFによる空間分解能の影響

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()

PSFによる中心部の強度の影響

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アニメーション作成コード

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の影響を意識することで、より正確なデータ解析が可能になります。本記事が、その理解の一助となれば幸いです。