あぷらなーと流・「非冷却」CMOSカメラのノイズ解析ごっこ

QHY5III174Mの問題

前に記事に書いた、人生初のモノクロカメラ QHY5III174Mですが、やはりというか懸念していた問題が浮上してきました。

・・・露光を伸ばすとアンプグローが出る・・・

さらにやっかいなことに、このアンプグローは非常に再現性がなく、撮影直後にそのままの状態で撮影したダークですら取り除くことができません。(つまり、温度だけでなく、そのほかの原因によっても変化している可能性があります)

一度こうなってしまうと、後処理で完全にアンプグローの影響を取り除くことはほぼ不可能です。そこでまずは、アンプグローの再現性がない理由を探るべく、このカメラのダークノイズ特性を調査しました。

ノイズ解析ごっこ

日本の天文界隈で大活躍され、あの天文雑誌『星ナビ』でも連載を書かれている、あのあぷらなーとさんが、『冷却CMOSカメラのノイズ解析ごっこ』と題して、ダークフレームの時間ノイズと空間ノイズを弁別する手法を検討されています。

この方法は、複数枚のダークフレームについて、各ピクセルの輝度値より、その中央値と標準偏差を計算するものです。
Fig.4は輝度中央値に対して標準偏差をプロットしたものです。図中に解説を書いてくださっていますので、詳しい説明は割愛しますが、まずは私のカメラで同じプロットを作成してみたいと思います。

問題

あぷらなーとさんは、プログラムの作成にMATLABという数値解析ソフトウェアを使用されています。しかし、このソフトは15,500円(個人ライセンス)もする高級ソフトです。プログラミングのまともな経験がない私が使いこなせるとも限らず、導入を断念しました。
そこで、MATLABの代わりに、無料で開発環境を構築できるプログラミング言語として、Pythonを選択しました。学生時代に少しだけ触れる機会があったこと、科学技術計算用のパッケージが豊富なことから選びました。
早速、プログラミングの勉強を兼ねて、あぷらなーとさんと同様の処理を行い、散布図を作成するプログラムの作成に取り組みました。

実装

今回作成したプログラムのソースコードを示します。(勉強を始めて一週間のため、よろしくない処理が入っているかもしれません。こうしたほうがいいよ！というのがありましたら、ぜひコメントいただけますと幸いです。)
興味のある方は、以下のコードをテキストエディタにコピーし、拡張子を.pyに変更し、Pythonがインストールされた環境で実行してみてください。

センサごとにピクセル数が異なるため、自動的にピクセル数 (=配列の行・列数)を読みに行く設定にしています。今回あぷらなーとさんにお手伝いいただいた際も、ピクセル数でハマったようですし・・・
(ここの処理はもっとスマートな書き方ができるだろうと思っています・・・)

"""

===========================================
CMOS camera noise analysis program ver.1.0
===========================================

*By: so-nano-car inspired by Mr. "Apranat"*

*License: BSD*

"""

import matplotlib.pyplot as plt
from matplotlib.pyplot import axis
import numpy as np
import pandas as pd
import glob
import astropy
import numpy as np
from astropy.visualization import astropy_mpl_style
from numpy.core.fromnumeric import shape, std
plt.style.use(astropy_mpl_style)
from astropy.utils.data import get_pkg_data_filename
from astropy.io import fits
import time

# Load start time
t1 = time.time()

# Load sensor pixel size (array size)
for img in glob.glob ('./hoge/*.fits'):
    imdata = fits.getdata(img, ext = 0)
array_y = np.shape(imdata)[0]
array_x = np.shape(imdata)[1]

# Load files and number of files: n, then store to 3-D numpy array
stack = np.empty((0,array_y,array_x))
for img in glob.glob ('./hoge/*.fits'):
    imdata = fits.getdata(img, ext = 0)
    stack = np.append(stack,imdata[np.newaxis,:],axis = 0)
count = np.shape(stack)[0]
print('{0} files loaded.'.format(count))

# Calculate median and std.dev. of each pixel
median = np.median(stack, axis = 0)
stddev = np.std(stack, axis = 0,ddof = 0)

# Reshape median and std.dev. array for plotting
x = median.reshape([array_y * array_x,1])
y = stddev.reshape([array_y * array_x,1])

# Export to csv file
data = np.concatenate([x,y], 1)
np.savetxt('./fuga/result.csv',data,delimiter = ',')

# Calculate elapsed time
t2 = time.time()
elapsed_time = t2 - t1
print("Elapsed time ={0} s.".format(elapsed_time))

# Plot results
plt.scatter(x,y,s = 0.3, marker = '.', color = "blue")
plt.grid(which = "both", linewidth = 0.5, alpha = 0.1)
plt.suptitle("Bias")
plt.xlabel('Median')
plt.ylabel('Std.dev.')
plt.xscale("log")
plt.yscale("log")
plt.xlim([100,100000])
plt.ylim([10,100000])
plt.show()

評価

解析対象となるダークフレームは、Fig.5 のように、できるだけ実際の撮影環境に近い形で撮影を行いました。ただし、室温は25℃前後あり、ライトフレームの撮影時の温度(7℃前後)より高く、熱ノイズは多いはずです。

結果

早速、作成した散布図を見てみましょう。ここは一つ邪道流に・・・

ででん！！

次に、Fig.4に倣って、自分なりの考察を入れてみました。

次に、番号付けした群に沿って、特徴を見ていきましょう。

①　オフセット

プロットの一番左端、輝度が低く、かつ標準偏差の値が小さい(ばらつきが小さい)ピクセルは、オフセット(バイアス)の存在を示唆していると考えられます。このことは、バイアスフレーム(10 us露光)の解析結果であるFig.8 からも読み取れます。

②　ダークノイズ群

プロットの中で一番特徴的な、直線的に伸びている群です。これは一般的に「ダークノイズ」と呼ばれている、熱電子による暗電流ノイズを拾ったピクセルであると考えられます。熱による電子の励起確率はポアソン分布に従うため、中央値(本来は平均値)と分散が等しい、つまり標準偏差は中央値の平方根にほぼ比例すると考えられます。(=両対数プロット上ではプロットが傾き0.5の直線となる)
実際にプロットから直線の傾きを概算してみると約0.57となり、これらの群はポアソン分布に従っていると言えそうです。(正常なダークノイズの挙動と言えそうです。)

③　放射線ヒット群

もっともピクセルが多い輝度値1000-2000付近の山から垂直に伸びる群です。これは放射線によるピクセルの感光であると考えられます。
あぷらなーとさんが、実際にこの群にいるピクセルを抽出し、時系列で輝度の変動を調査してくれました。(大変そうだったのでお願いしていなかったのですが・・・本当にありがとうございました。)

Fig.9より、例えば[1098,220]という座標のピクセルでは輝度値が1フレームだけ異常に大きい値を吐き出しています。このことから、この[1098,220]ピクセルは当該フレームを撮ったときだけ宇宙線による影響を受けたと考えられます。
一方で、[165,1166]ピクセルでは宇宙線ヒットに類似した輝度変動だけでなく、ランダムな輝度の変動が見られます。この原因については明らかになっておらず、今後調査予定です。

④　酩酊ピクセル

②群の上側に点在するピクセルは、熱ノイズが従うはずのポアソン分布に「従わない」ピクセル、つまり異常ピクセルであると考えられます。Fig.9でも言及されていますが、輝度がランダムに変動するピクセルに対して、あぷらなーとさんが「酩酊ピクセル」と命名したものです。

⑤　謎の分裂した群

②の上に、直線状に伸びる分裂した群が見られます。②でも述べたように、「正常な」ダークノイズは、傾きが0.5の直線上に乗るはずです。(Fig.10中の赤点線)
しかし、この群⑤は傾きが大きくなっています。そこで、あぷらなーとさんが、撮影開始直後から64コマ、撮影開始10分後から64コマに分けて解析をしてくださいました。

Fig.11 より、群⑤に位置するピクセルは、その輝度が時間経過とともに低くなり、安定していることが分かります。一方、群②に位置するピクセルは、その輝度が時間経過にかかわらずほぼ一定となっています。
このことから、センサ内のピクセルが受ける熱ノイズの影響はピクセル間で一定でなく、ピクセルごとに異なっている可能性が示唆されました。この原因については、今後の課題ですね。

結論

今回の記事では、QHY5III174Mのダークフレームについて時系列解析を行い、ダークフレームの特性をつかむことを試みました。結果として、温度による影響がセンサ内でも一様でなく、このことがアンプグローの再現性低下に寄与している可能性があります。

このことからいえるのは、このカメラで撮影をする際は、アンプグローの影響を極力受けない露光時間以下で、多数枚の露光を重ねるのがよさそう、ということだと思います。

蛇足

Fig.6に示したグラフは、センサごとに大きく異なる形を示します。Twitter上のタイムラインでは、このプロットが「ナイキ」のマークに見えると話題でした。みなさんのカメラも、同様の解析をすると、面白い模様が見えてくるかもしれません。

最後に、解析の原案を作成され、私のデータを使って解析までしていただいたあぷらなーとさんに、重ねてお礼申し上げます。ありがとうございました。