本稿では、Linux系OSのシェルとIRAFを用いて「淡々と」、撮像観測データを画像処理する方法を述べる。勿論、もっとエレガントに行うことも可能だが、ここではあくまでも初心者向きに「泥臭く」「地味に」行うことをモットーとする。

対象とする観測画像は、東京大学木曽観測所の105cmシュミット望遠鏡と2kCCDカメラで得られた撮像データである。2kCCDカメラは既に2011年に次世代観測装置（KWFC）に取って変わられているが、観測で得られたデータは現在でも SMOKA などを通して得ることが可能である。

2kCCDカメラで得られたデータは、「kcd」のあとに 5 ないし 6 ケタの通し番号を付けたFITS形式画像で保存される。つまり、ファイル名は「kcdxxxxxx.fits」（xxxxxxは通し番号）となる。リダクションをはじめとする観測データの画像処理では、作業を進めるたびに新たな画像データが作成されるが、ここでは、通し番号を維持したまま、ファイル名冒頭の「kcd」を他の記号と入れ替えて新しい画像ファイルを作成する。また、作成される画像ファイルの全ては、生データと同じ場所に置くものとする。もしも、画像ファイルを特定のディレクトリに保存したい場合などは各自で工夫して欲しい。

なお、本稿では Linux のシェルのプロンプトを「$>」、IRAFのプロンプトを「cl>」と表記することにする。

• 第１部 --- １次処理（整約：リダクション）
  • 01. 簡易観測ログの作成
  • 02. バイアス画像の合成
  • 03. バイアス成分のレベル合わせと除去
  • 04. フラット画像の規格化
  • 05. フラット画像の作成
  • 06. 天体画像のフラット・フィールディング
  
  
• 第２部 --- ２次処理
  • 07. バッド・ピクセルの補正
  • 08. オーバースキャン領域（ OSR = Overscan Region ）の除去
  • 09. 背景光（スカイ）成分の除去
  
  
• 第３部 --- 天体画像の合成 冒頭部能書きのみ
  • 10. 星像サイズとフラックス・レベルの測定 作成中
  • 11. 星像サイズの調整 作成中
  • 12. フラックス・レベルの調整 作成中
  • 13. 画像間の位置合わせ 作成中
  • 14. 天体画像の合成 作成中

観測で得られた直後で、まだどのような処理も施されていないデータを生データ（raw data）と呼ぶ。生データの信号 I_raw の中には、天体からの信号 I_obj以外にも、次のような様々な信号が含まれている。初めに 1) 星野光と呼ばれる空間的に分解できない無数の遠方天体による背景光、次に、2) 黄道光と呼ばれる太陽系円盤に存在するチリによって散乱された太陽光による背景光、続いて、3) 大気光と呼ばれる高層大気の発光による背景光、の３成分でこれら３者はまとめて「夜天光」と呼ばれる。背景光には 4) 大気による都市光（city light）の反射光も含まれる。夜天光と都市光は併せて、単に「背景光（background light）」または「スカイ（sky）」I_skyと呼ばれる。これら 1) から 4) は宇宙・地球・観測環境に起因する信号であるが、さらに、観測装置自身に起因する信号も存在する。まず、5) 光子ではなく、熱による暗電流信号（dark current：ダーク）I_dark。3) 出力の雑音情報を落とさないために、人為的に付加されるバイアス信号（bias）I_bias、である。そして、信号そのものでは無いが、CCDチップの各ピクセル（画素のこと）間の感度は必ずしも同じではないため、これによって感度むら Sが生じる。これら生データに含まれる信号の情報をまとめると、

I_raw = （ I_obj ＋ I_sky ）× S ＋ I_bias ＋ I_darkとなる。

従って、観測によって得られた生データから、天体由来の信号I_objを取り出すためには暗電流（ダーク）成分I_darkとバイアス成分I_biasを取り除き、感度ムラSを補正した後、背景光成分I_skyを除去すれば良いことになる。一般に、感度ムラ補正までを１次処理（整約：リダクション）と呼ぶ。

なお、共同利用に供される可視光CCDカメラの多くでは、暗電流成分I_darkは無視できるほど小さいことが多い。2kCCDにおいてもI_darkの処理は省略することが普通である。

各FITSファイルのヘッダーを参照することで、簡易な観測ログを作成する。このログには、後々画像解析に必要となる画像の典型的な輝度の値とオーバースキャン領域（OSR）の値も記載する。

まず、シェルと IRAF のカレント・ディレクトリを、生データが置いてあるディレクトリにする。そこで以下の一連の作業を実行する。

$> ls -1 kcd*.fits > ___raw.lst [RET] ： 全生データのリスト作成 cl> hselect @___raw.lst FRAME-ID,DATA-TYP,OBJECT,FILTER01,EXPTIME,JST,ZD yes > ___rawhslct.tmp [RET] ： FITSヘッダーの読み取りと出力 $> awk '{{if(substr($3,1,3)=="kcd"){print $3" "$4" "$5" "$6" "$7" "$8" "$9}}{if(substr($1,1,3)=="kcd"){print $0}}}' ___rawhslct.tmp > ___rawhslct.dat [RET] ： FITSヘッダー出力の整形と出力 $> awk '{print $1"[11:2038,1:2048]"}' ___raw.lst > ___raw0.lst [RET] ： 非オーバースキャン領域のリスト作成 cl> imstat @___raw0.lst fields=midpt format=no > ___rawmidpt.tmp [RET] ： 非オーバースキャン領域の輝度の代表値(メジアン値)出力 $> awk '{print $1"[3:8,3:2046]"}' ___raw.lst > ___rawosr.lst [RET] ： オーバースキャン領域のリスト作成 cl> imstat @___rawosr.lst fields=midpt format=no > ___rawosr.tmp [RET] ： オーバースキャン領域の輝度の代表値(メジアン値)出力 $> paste ___rawhslct.dat ___rawmidpt.tmp ___rawosr.tmp > raw.obslog [RET] ： FITSヘッダー情報、非OSR情報、OSR情報のマージと出力 $> rm ___*.* [RET] ： 不要なファイルの削除

これで、１行に、FITSファイル名、データタイプ、観測対象名、使用フィルター名、露光時間、観測開始時刻、観測開始天頂距離、典型的な輝度の値、OSRの代表値、が記された簡易観測ログが、「raw.obslog」という名前で出力される。このファイルの情報を元に、以下の画像処理を進めることになる。

前節で作成した簡易観測ログ「raw.obslog」を元に、バイアス画像を選出し、それを合成することでより高精度のバイアス画像を作成する。

$> awk '{if($2=="BIAS" && $8<4350){print $1}}' raw.obslog > ___rawbias.lst [RET] ： 平均の輝度が4350(ADU)未満のバイアス画像のリストを出力 cl> imcombine @___rawbias.lst BIAS.fits combine=average reject=none offsets=none [RET] ： バイアス画像の合成（平均） cl> imdel @___rawbias.lst [RET] ： 不要なファイルの削除。ただし、生データを保存する場合、この作業はスキップ。 $> rm ___rawbias.lst [RET] ： 不要なファイルの削除

これで平均的なバイアス画像が、「BIAS.fits」というファイルで作成される。

前節で作成したバイアス画像「BIAS.fits」を用いて、フラット画像と天体画像からバイアス成分を取り除く。この際に、OSRの値を考慮して、より高精度でのバイアス成分除去を行う。なお、2kCCDでは、FITSヘッダー「DATA-TYP」によるドームフラットとスカイフラット、トワイライトフラットなどの区別はされておらず、一律に「FLAT」と記載されている。

cl> imstat BIAS.fits[3:8,3:2046] fields=midpt format=no > biasosr.dat ： 合成バイアス画像のOSRの値を出力 $> awk '{if(($2=="FLAT" || $2=="OBJECT") && $8<40000){print $1" "$8" "$9}}' raw.obslog > ___rawflatobj.log [RET] ： 平均の輝度が40000(ADU)未満のフラット画像/天体画像のファイル名と輝度の値、OSRの値を出力 $> awk '{if(file==1){biasosr=$1;};if(file==2){print $1" "$3/biasosr}}' file=1 biasosr.dat file=2 ___rawflatobj.log > ___rawflatobj.tmp [RET] ： フラット画像名/天体画像名と"これら画像のOSRの値を合成バイアス画像のOSRの値で規格化した"値を出力 $> awk '{print $1}' ___rawflatobj.tmp > ___rawflatobj.lst ： フラット画像/天体画像の生データのリスト出力 $> awk '{print $2}' ___rawflatobj.tmp > ___flatobj2bias.dat ： レベル合せ用定数のリスト出力 $> awk '{print "bias"substr($1,index($1,"kcd")+3)}' ___rawflatobj.tmp > ___flatobj2bias.lst ： レベル合せ済みバイアス画像名のリスト出力 cl> imarith BIAS.fits * @___flatobj2bias.dat @___flatobj2bias.lst [RET] ： レベル合せ済みバイアス画像の作成 $> awk '{print "bs"substr($1,index($1,"kcd")+3)}' ___rawflatobj.tmp > ___bsflatobj.lst ： バイアス差し引き済みフラット画像/天体画像のリスト出力 cl> imarith @___rawflatobj.lst - @___flatobj2bias.lst @___bsflatobj.lst [RET] ： 生フラット画像/生天体画像からレベル合せ済バイアス画像を差し引く cl> imdel @___rawflatobj.lst [RET] ： 不要なファイルの削除。ただし、生データを保存する場合、この作業はスキップ。 cl> imdel @___flatobj2bias.lst [RET] ： 不要なファイルの削除。ただし、中間データを保存する場合、この作業はスキップ。 $> rm ___*.* [RET] ： 不要なファイルの削除 $> rm biasosr.dat [RET] ： 不要なファイルの削除 $> rm BIAS.fits [RET] ： 不要なファイルの削除

これで、OSRを考慮した上でバイアス成分を取り除いたドームフラット画像と天体画像が、「bsxxxxx.fits」（xxxxxは通し番号）というファイル名で作成される。

簡易観測ログ「raw.obslog」を元に、バイアス成分除去後のフラット画像を選出し、それらを非OSRの輝度の代表値で除することで規格化する。これによって全てのフラット画像のレベルを「平均１」に合わせる。

$> awk '{if($2=="FLAT"){print "bs"substr($1,index($1,"kcd")+3)}}' raw.obslog > ___bsflat.lst [RET] ： バイアス成分除去済みのフラット画像のリストを出力 $> awk '{if($2=="FLAT"){print "bs"substr($1,index($1,"kcd")+3)"[11:2038,1:2048]"}}' raw.obslog > ___bsflatnoosr.lst [RET] ： バイアス成分除去済みのフラット画像の非OSRリストを出力 cl> imstat @___bsflatnoosr.lst fields=midpt format=no > ___bsflatnoosr.dat [RET] ： バイアス除去済みフラット画像の輝度の代表値の出力 $> awk '{if($2=="FLAT"){print "nm"substr($1,index($1,"kcd")+3)}}' raw.obslog > ___nmflat.lst [RET] ： 規格化したフラット画像のリストを出力 cl> imarith @___bsflat.lst / @___bsflatnoosr.dat @___nmflat.lst [RET] ： バイアス除去済みフラット画像の規格化 cl> imdel @___bsflat.lst [RET] ： 不要なファイルの削除。ただし、中間データを保存する場合、この作業はスキップ。 $> rm ___*.* [RET] ： 不要なファイルの削除

これでバイアス成分を取り除き、レベル合せが完了したドームフラット画像が、「nmxxxxx.fits」（xxxxxは通し番号）というファイル名で作成される。

フィルターバンド毎にフラット画像を合成する。簡易観測ログ「raw.obslog」を元に、規格化されたフラット画像をフィルター毎に選出し、メジアンフィルターを通して合成する。ここでは V バンドを取り上げて説明するが、2kCCDでは他にも広帯域撮像用として「U」「B」「R」「I」フィルターが、狭帯域撮像用として「N487」「N499」「N519」「Ha6417」「Ha6577」「Ha6737」フィルターなどが用意されていた。なお、基本的に、2kCCD の U, B, V は Johnson システム、R, I は Cousisns システムに準拠したものとされるが、実際には U および I バンドのフィルター関数は、推奨されているものからのズレが大きめである。

$> awk '{if($2=="FLAT" && $4=="V"){print "nm"substr($1,index($1,"kcd")+3)}}' raw.obslog > ___nmflatV.lst [RET] ： 規格化した V バンドのフラット画像のリストを出力 cl> imcombine @___nmflatV.lst DFlatV.fits combine=median reject=none offsets=none [RET] ： V バンドフラット画像の合成 cl> imdel @___nmflatV.lst [RET] ： 不要なファイルの削除。ただし、中間データを保存する場合、この作業はスキップ。 $> rm ___*.* [RET] ： 不要なファイルの削除

この作業で、V バンドのフラット画像が「DFlatV.fits」という名前で作成される。明記していなが、ファイル名冒頭の「D」はドームフラット（Domeflat）であること、「V」は V バンドを示している。

ピクセル間の感度ムラ補正を「フラット・フィールディング」という。フィルターバンド毎に、バイアス除去後の天体画像のフラット・フィールディングを行う。簡易観測ログ「raw.obslog」を元に、バイアス除去後の天体画像を選出し、対応するフラット画像を除することでフラット・フィールディングを行う。ここでも、前節に引き続き V バンドを取り上げて説明する。

$> awk '{if($2=="OBJECT" && $4=="V"){print "bs"substr($1,index($1,"kcd")+3)}}' raw.obslog > ___bsobjV.lst [RET] ： バイアスを除去した V バンドの天体画像のリストを出力 $> awk '{if($2=="OBJECT" && $4=="V"){print "ff"substr($1,index($1,"kcd")+3)}}' raw.obslog > ___ffobjV.lst [RET] ： フラット・フィールディング後の V バンドの天体画像のリストを出力 cl> imarith @___bsobjV.lst / DFlatV.fits @___ffobjV.lst [RET] ： V バンド天体画像のフラット・フィールディング cl> imdel @___bsobjV.lst [RET] ： 不要なファイルの削除。ただし、中間データを保存する場合、この作業はスキップ。 $> rm ___*.* [RET] ： 不要なファイルの削除

これで、フラット・フィールデイングが完了し、感度ムラ補正が行われた画像が「ffxxxxx.fits」（xxxxxは通し番号）という名前で作成される。これで、１次処理（整約：リダクション）は終了である。

フラット・フィールディングまでの処理を１次処理や整約、リダクションと呼ぶのに対して、フラット・フィールディング以降の画像処理を「２次処理」と言うらしいのだが、正直、本稿の筆者はそのような表現を聞いたことが無い。若い方々の大部分は、「観測データの画像処理全般をリダクションと呼ぶ」と思っているのでは無いだろうか？（笑）この第２部では、１枚の撮像画像をサイエンス・データとして完成させるところまでを行うことにする。

手順としては、フラット・フィールディングの後に背景光の除去となるが、その前にここでは、バッド・ピクセルの補正を行う。バッド・ピクセルとはフラット・フィールデングでは補正し切れないような欠陥ピクセルのことである。常に極端に暗く、または、極端に明るく応答し、これらが直線上に配列していることも多く、画像中で明るい/暗い線として現れる。バッド・ピクセルが線状に現れたものはバッド・カラムと呼ばれることもある。基本的に画像中で、これらバッド・ピクセル、バッド・カラムの位置が変化するとは無いはずだが、経年劣化によって新たなバッド・ピクセルやバッド・カラムが現れることはあるようである。

なお、バッド・ピクセルの補正に用いる IRAF/fixpix は、画像を上書きしてしまうので、事前に補正前の画像のバックアップを取るなど注意して使用して欲しい。

$> ls -1 ff*.fits > ___ffobj.lst [RET] ： フラット・フィールディング済みの天体画像のリスト出力 $> vi badpix2k.dat [RET] ： バッド・ピクセルの位置のリスト作成（手作業による） 以下は 2kCCD のバッド・ピクセル補正用の位置情報（copy & paste して利用して欲しい） 669 669 659 769 377 377 803 828 703 703 824 834 712 712 131 232 1022 1022 191 210 1358 1360 1 2048 940 940 754 767 969 969 716 726 837 839 920 2048 1180 1180 973 994 1204 1204 918 1205 1233 1233 1093 1132 1325 1325 994 1155 1314 1314 1012 1032 1335 1335 974 984 1386 1386 1016 1055 803 804 1976 2048 1097 1098 2034 2048 1649 1649 2023 2048 1705 1705 1763 1788 1993 1993 1667 1752 1655 1655 1590 1636 1423 1423 1167 1202 2022 2022 724 741 1486 1486 601 624 1564 1564 199 208 883 883 1204 1242 883 883 1307 1420 1204 1204 1366 1537 1325 1325 1500 1750 1285 1285 1048 1240 917 917 1031 1069 883 883 495 590 892 892 608 623 969 969 726 749 1012 1012 191 207 1303 1303 1027 1062 1651 1651 1966 1998 1993 1993 1752 1872 cl> fixpix @___ffobj.lst masks=badpix2k.dat verbose=no [RET] ： 周辺ピクセルからの内挿によるバッド・ピクセルの補正 $> rm badpix2k.dat [RET] ： 不要なファイルの削除。ただし、バッド・ピクセルの位置データを保存する場合、この作業はスキップ。

背景光除去に先立って、画像データの両端に付加されたオーバースキャン領域（ = OSR ）を除去する。正確には、天体が撮影された非OSRを別のファイル名でコピーすることになる。

$> ls -1 ff*.fits > ___ffobj.lst [RET] ： フラット・フィールディング済み、かつ、バッド・ピクセル補正済みの天体画像のリスト出力 $> awk '{print $1"[11:2038,1:2048]"}' ___ffobj.lst > ___ffobjnoosr.lst [RET] ： 天体画像の非OSRのリスト出力 $> awk '{print "oc"substr($1,index($1,"ff")+2)}' ___ffobj.lst > ___ocobj.lst [RET] ： OSR除去後の天体画像のリスト出力 cl> imcopy @___ffobjnoosr.lst @___ocobj.lst [RET] ： OSRの除去 cl> imdel @___ffobj.lst [RET] ： 不要なファイルの削除。ただし、中間データを保存する場合、この作業はスキップ。 $> rm ___*.* [RET] ： 不要なファイルの削除

これで、オーバースキャン領域を除去された天体画像が「ocxxxxx.fits」（xxxxxは通し番号）という名前で作成される。

フラット・フィールディング、バッド・ピクセルの補正、OSRの除去が完了した天体画像から、背景光（スカイ）成分を、多項式の２次元面でフィッティングすることで取り除く。

$> ls -1 oc*.fits > ___ocobj.lst [RET] ： 背景光除去前の天体画像のリスト出力 $> awk '{print "ss"substr($1,index($1,"oc")+2)}' ___ocobj.lst > ___ssobj.lst [RET] ： 背景光除去後の天体画像のリスト出力 cl> imsurfit @___ocobj.lst @___ssobj.lst xorder=2 yorder=2 type_ou=residual functio=leg cross_t=no xmedian=64 ymedian=64 [RET] ： 天体画像からの背景光除去。平面フィットにはクロスタームが無い（xy項の係数が０）２次式を用い、補間にはルジャドル、フィッティングに用いる領域は64x64である。 cl> imdel @___ocobj.lst [RET] ： 不要なファイルの削除。ただし、中間データを保存する場合、この作業はスキップ。 $> rm ___*.* [RET] ： 不要なファイルの削除

これで、背景光を除去された天体画像が「ssxxxxx.fits」（xxxxxは通し番号）という名前で作成される。これで１枚の撮像データとしては、画像処理が完了したことになり、「完成」したことになる。

一枚の撮像データに対する画像処理としては、第２部までの手続きで事実上、完了・完成している。しかし実際には、複数の撮像データを合成することで、より高精度な画像データの作成を行うケースも頻繁に存在する。ところがこの場合、撮像データ間で様々な条件が異なることが多いため、単純にそのまま合成する訳には行かない。考慮・調整すべき条件は、1) 主に大気の乱れに起因する星像サイズの違い、2) 大気の水蒸気量・透明度に由来するフラックス・レベル（輝度レベル）の違い、3) 望遠鏡の指向精度やディザリング観測による位置のズレ、の３者である。