電波干渉計の観測を雑にPythonで表現してみたという記事で電波干渉計の原理をPythonで非常に単純化した形で表現しました。この記事では続編として、ダーティイメージから「汚れ」のない画像を復元する手法について、引き続きPythonで表現していきます。

なお、ここでは添付の画像を生成するコードは示していません。添付画像の生成を含むPythonコードの全体像はGistにまとめました。

環境

以下のコードは下に示す環境で実行しました。

Python 3.10.8
numpy 1.23.4
scipy 1.9.3
scikit-learn 1.1.3

また、添付の画像はmatplotlibで生成しました。

電波干渉計のおさらい

前回記事のまとめより、電波干渉計とは以下のような観測装置です。

電波干渉計は天体画像のフーリエ成分を観測する「空間周波数フィルタ」
なるべく多くのフーリエ成分を観測することで、フーリエ成分から天体画像を復元する
ただしデータが欠損するので画像の復元には工夫が必要

ここでは簡単のため2次元画像の代わりに1次元配列を使ってアルゴリズムを実装します。なのであらためて、1次元配列として入力画像とそのフーリエ成分を生成しておきます。なお、ここでは画像のフーリエ成分を「ビジビリティ」と呼ぶことにします。電波干渉計はビジビリティを測定する装置ということになります。

import collections

import numpy as np

# 元画像（単純化のため1次元配列）
nx = 128
image_original = np.zeros(nx, dtype=float)

# 2つの放射源
LineProperty = collections.namedtuple('LineProperty', ['position', 'width', 'intensity'])
line_list = [
    LineProperty(position=60, width=2, intensity=10),
    LineProperty(position=70, width=4, intensity=3),
]
for line in line_list:
    line_start = line.position - line.width // 2
    line_end = line_start + line.width
    image_original[line_start:line_end] = line.intensity

# 電波画像のフーリエ成分をビジビリティと呼ぶ
visibility = np.fft.fft(image_original, norm='ortho')

ビジビリティを完全に測定することはできず、データは必ず欠損します。データの欠損状況を表すマスクを、UVカバレッジと呼びます。UVカバレッジを逆フーリエ変換すると、中心にピークを持ったビームのような分布が得られます。これが電波干渉計の点像分布関数（PSF）になります。PSFにはダーティイメージに現れる偽の構造の特徴が含まれています。その意味で電波干渉計のPSFをダーティビームと呼びます。ダーティイメージの「汚れ」を全てPSFに押しつけたという見方ができると思います。ダーティビームの主ビーム成分をガウシアン近似したものは、ダーティビームとの対比でしょうか、クリーンビームと呼びます。

import scipy.signal as signal

# フーリエ成分の観測が不完全であることを表すマスク
# これをUVカバレッジと呼ぶ
uv_limit = nx // 4
mask = np.ones_like(visibility, dtype=int)

# 観測は最大のアンテナ間隔で制限されている
mask[uv_limit:-uv_limit] = 0

# 原理的に観測可能な範囲でも観測は不完全
zeros = np.random.randint(1, uv_limit, size=uv_limit // 2)
mask[zeros] = 0
mask[zeros * -1] = 0
mask[1:3] = 1
mask[0] = 0
mask[-2:] = 1

# ダーティビーム = UVカバレッジの逆フーリエ変換
# 干渉計のPSFに相当する
psf = np.fft.fftshift(np.fft.ifft(mask, norm='ortho').real)
psf /= psf.max()

# クリーンビーム
# 本来はダーティビームの主ビーム成分をガウシアンでフィットする
clean_beam = signal.windows.gaussian(nx, std=1.2)
clean_beam /= clean_beam.max()

不完全なビジビリティ観測の欠損部分をゼロ埋め（ゼロパディング）するとダーティイメージが得られるのでした。

visibility_observed = visibility * mask
visibility_observed_shift = np.fft.fftshift(visibility_observed)

# ダーティイメージ
# 観測されなかったビジビリティをゼロパディングして逆フーリエ変換することで得られる
image_observed = np.fft.ifft(visibility_observed, norm='ortho').real

電波干渉計の画像復元法その１: 伝統と信頼のCLEAN法

ダーティイメージの「汚れ」をPSFに押しつけることにすると、ダーティイメージは真の画像に汚れたPSFであるダーティビームが畳み込まれた（convolve）結果と捉えることができます。したがって、ダーティイメージから汚れたPSFを取り除く（deconvolve）ことができれば真の画像が得られると期待できます。このような考え方から、伝統的にダーティイメージから真の画像を得る手法をデコンボリューション（deconvolution）と呼んでいます（と、自分なりに勝手に解釈しています）。デコンボリューションの代表的なアルゴリズムはCLEAN（クリーン）法です。汚れを取るからCLEANなのでしょうか。

CLEAN法の原理を以下にPythonで示します。プロが見たら怒られそうなくらい単純化しているのですが、本質は捉えていると信じています。残差画像、クリーン成分、モデル画像の３つを準備します。手順は以下の通りです。

残差画像にダーティイメージををセットする。クリーン成分とモデル画像は全てのピクセルにゼロをセットする。
残差画像から絶対値が最大のピクセルを探す
そのピクセルの値にゲイン（通常10%程度）をかけた値をクリーン成分に登録
クリーン成分にクリーンビームをたたみ込んだ画像をモデル画像とする
クリーン成分にダーティビームをたたみ込んだ画像を残差画像から差し引く
上記2〜5を収束するまで繰り返す

この手順は、一言でいえばダーティビームのクリーンビームへの置き換えです（手順4と5）。今、画像の「汚れ」は全てダーティビームに押し付けているので、ダーティビームをクリーンビームに置き換えてやれば「汚れ」のないクリーンな画像が得られるだろうという考え方と理解しています。実際この手順を繰り返すことで残差画像（＝ダーティイメージ）の汚れが取れた画像がモデル画像として得られていることがわかります。

import scipy.signal as signal

# ごく単純なCLEAN法の実装

# ゲインファクター
# 数値的な振動を抑制するために、ちょっとずつ処理を進める
alpha = 0.1

# CLEAN成分を保持する配列
clean_component = np.zeros_like(image_observed)

# モデル画像 = CLEAN成分をクリーンビームでたたみ込んだ画像
model_image = signal.convolve(clean_component, clean_beam, mode='same')

# 残差画像
residual = image_observed - signal.convolve(clean_component, psf, mode='same')

# CLEAN法の実装
for i in range(10):
    # ピーク位置を探す
    peak_position = np.argmax(residual)

    # CLEAN成分に見つけたピークをゲインファクターをかけて登録する
    clean_component[peak_position] += residual[peak_position] * alpha

    # モデル画像の更新
    model_image = signal.convolve(clean_component, clean_beam, mode='same')

    # 残差画像の更新
    # CLEAN成分にダーティビームをたたみ込んだ画像を差し引く
    residual = image_observed - signal.convolve(clean_component, psf, mode='same')

ここでは上から順に繰り返し回数0回（初期状態）、10回、50回の結果を示しました。繰り返し回数を増やすにしたがって、モデル画像（右列）において入力画像の２つの放射成分が再現され、しかも偽の構造がほぼ消えていることがわかります。CLEAN法では、モデル画像に残差画像を加えることで最終結果とするのが慣例です。

# CLEANの最終画像 = モデル画像 + 残差画像
clean_image = model_image + residual

これで、不完全なフーリエ成分の観測から入力画像に近い結果が得られました。強度の問題はとりあえず無視すると、見事に２つの放射成分を再現できました。

電波干渉計の画像復元法その2: 新興勢力RML法

CLEAN法は実用面でのさまざまな改良が現在でも続いており、電波干渉計業界ではほぼデファクトスタンダードです。いっぽうで、最近電波干渉計の画像復元においてもデータ科学的な手法の開発が非常にホットな話題になっています。電波干渉計の画像復元は、欠損したデータから画像を復元するという意味で劣決定問題になっています。したがって、データ科学的手法（統計学的手法）によって真の画像を復元するというアプローチも可能です。この方向性での画像生成手法、特に制限付き最尤推定法（RML法）による画像生成の手法は、ブラックホールシャドウの撮影に成功したイベントホライゾンテレスコープで脚光を浴びました。

RML法では、事前知識にもとづく制約を導入することにより、解けない問題を解けるようにする手法です。ここではRML法の一例として、L1ノルム（絶対値の和）の制約を課して問題を解くLASSO回帰によって画像を復元します。LASSOの制約はL0ノルム（ゼロでない要素の個数）の代替としてL1ノルムを使っているもので、制約を強めれば強めるほどゼロの多い解（スパースな解）が得られるような仕組みになっています。

from sklearn import linear_model

# 制限付き最尤推定法（RML）の実装
# ここではLASSOにより画像を復元する

coeff = []
vis2 = []

# 観測されたフーリエ成分だけを処理対象とする
# ゼロパディング不要であることに注意
components = np.arange(nx)[mask != 0]

# フーリエ変換係数を配列に保存する
for k in components:
    # 複素数版のフーリエ変換係数は以下のようにかける
    # 
    # _coeff = np.exp([complex(0, x) for x in range(nx)]) * -2 * np.pi * k / n
    #  
    # しかし、scikit-learnのLASSOは複素数配列をサポートしていないため、係数を実部と虚部に分ける
    # 
    # フーリエ変換係数の実部
    _coeff_real = np.cos(np.arange(nx) * 2 * np.pi * k / nx) / np.sqrt(nx)
    # フーリエ変換係数の虚部
    _coeff_imag = -np.sin(np.arange(nx) * 2 * np.pi * k / nx) / np.sqrt(nx)
    coeff.append(_coeff_real)
    coeff.append(_coeff_imag)
coeff = np.asarray(coeff)

# 観測されたフーリエ成分を抜き出す
observed = []
for k in components:
    observed.append(visibility[k].real)
    observed.append(visibility[k].imag)
observed = np.asarray(observed)

# LASSOソルバーのセットアップ（制約項の強度パラメータを0.01に設定）
clf = linear_model.Lasso(alpha=0.01)

# 解く
clf.fit(coeff, observed)

# 得られた画像を取得する
lasso_image = clf.coef_

CLEAN法と同じく２つの放射成分を見事に再現することができました。

CLEAN法とRML法の比較

いずれの手法を採用するにしても、もともとの問題が劣決定である以上、真の画像に辿り着くことはできないことには注意が必要です。あくまで「もっとも確からしい画像」を得ることが目的であり、結果の解釈は復元方法の性質を踏まえて慎重に行う必要があります。

CLEAN法の強みはやはり長年の知見の積み重ねによる安定性です。標準的な設定にのっとって画像復元する限りにおいては、処理方法にケチをつけられることはないでしょう。一方で、CLEAN法では繰り返し処理の打ち切り条件にやや不定性があります。また、原理的にクリーンビーム以上の解像度を持った画像が得られないという弱点があります。

一方RML法は、新しい手法だけに結果の解釈はClEAN法以上に慎重さが求められます。また、ハイパーパラメータの不定性があります。上記に示したLASSOの場合、ソルバーのパラメータとしてαを与えました。当然αの値によって結果は変わります。では数あるαの中で、どの値を取るべきでしょうか。何らかの方法で最適なハイパーパラメータを選ばなければなりません。このようにいくつか課題があるとはいえ、RML法ではCLEAN法に無い強みがあります。ひとつはクリーンビームのような解像度の制限がないことです。原理的には無限の解像度を達成できますが、科学的に意味があるのはクリーンビームの数倍程度と言われています。ですが仮にCLEAN法に比べて2倍の解像度で画像が得られたのだとしても、「たかが2倍、されど2倍」で科学的には大きな成果です。もうひとつの強みは、ゼロパディングの影響がないことです。CLEAN法は画像ドメインで処理をおこなうため、どうしてもビジビリティを逆フーリエ変換しなければならず、ゼロパディングの影響は不可避でした。一方RML法はビジビリティドメインで問題を解くため、ゼロパディングの影響は受けません。これは大きな強みです。

ALMA望遠鏡のための画像復元ソフトウエアPRIISM

私たちはRML法の流れに乗ってALMA望遠鏡のための画像復元ソフトウエアPRIISMを開発しました。LASSOではL1ノルムを制約項としていましたが、この制約だけでは常にスパースな解が選ばれてしまいます。ALMAではさまざまな天体が観測対象であり、中にはぼんやりと雲のように広がった天体もあります。そのような天体に対しても適用可能なソフトウエアにするために、PRIISMではL1ノルムに加えてTotal Squared Variation (TSV)という制約を追加しています。TSVは隣り合うピクセル間の値の差の二乗和で、TSVによる制約を強めると「滑らかな解」が得られる傾向になります。PRIISMではスパースな解を好むL1制約と滑らかな解を好むTSV制約の競合によって最適な解が決まります。このとき最適な制約の強さ（ハイパーパラメータ）は交差検証という方法で決定します。交差検証は機械学習の汎化性能を測定する際などにも用いられる手法で、得られた解が観測データに過剰適合しないようにチェックするという意味合いを持たせています。PRIISMはまだ発展途上ではありますが、PRIISMを使って論文を書いてくれる研究者も出てきました。CLEAN法を超えるとまではいきませんが、この業界で一定の地位が得られればいいなと思いながら細々と開発を続けています。

以上、まとめます。

電波干渉計の画像復元は、CLEAN法がデファクトスタンダードである
最近制限付き最尤推定法（RML法）により画像を復元する手法が流行っている
筆者らはRML法の流れをくむPRIISMというソフトウエアを開発した