すばる中秋の名月の学校 2015
模擬研究 太陽系コース
担当： 寺居 剛（国立天文台ハワイ観測所）
2015 年 9 月 28 日 – 10 月 1 日
1
はじめに
太陽系小天体（Small Solar System Bodies）とは惑星・準惑星・衛星を除いた太陽を公転する天体の総
称で、小惑星や彗星などが該当する。中でも、火星軌道と木星軌道の間の小惑星帯にはメインベルト小
惑星（main-belt asteroids）が、海王星軌道以遠の領域には太陽系外縁天体（trans-neptunian objects；
TNOs）が群集している。
太陽系小天体は地球との相対速度によって天球上を時々刻々と動いて見える。そのため、太陽系小天体
を見つけるには、同じ領域を時間をおいて複数回撮像し、背景の恒星などに対して動いている天体（移
動天体）を探せばよい。
天球上での移動速度は天体の軌道および地球との位置関係によって異なる。例えば衝（太陽の反対点）
の位置を観測した場合、典型的な移動速度はメインベルト小惑星で ∼ 30–45 arcsec/hr、外縁天体で ∼
3 arcsec/hr である。逆に言えば、移動速度を測定することで天体の大まかな軌道を推定することができ
る。ただし、衝から離れるほど軌道の違いによる速度差が小さくなり、分離は困難になる。
多くの太陽系小天体は黄道面（ecliptic plane）に沿った軌道を持つため、移動方向は黄道の向きに近い
ものが多い。また、一般的には地球の方が公転速度が大きいため、観測領域が衝から大きく（メインベ
ルト小惑星なら 50◦ 程度）離れていなければ、小天体は負の黄経方向に動く。
すばる望遠鏡＋ Suprime-Cam による観測では、高い集光力と広い視野を生かして多くの太陽系小天体
を検出することが可能である。黄道面領域を撮像すると、１視野で 100 個程度の小惑星が検出される。
それらのほとんどは未知の天体である。
1
使用するデータ
2
2.1
ターゲット領域
本コースでは以下のデータを解析する。
EXP ID： SUPE01354310 – SUPE01354350
観測日： 2012 年 7 月 25 日（UT）
対象天体： RXCJ2034.9-2143 (赤方偏移 z = 0.19 の銀河団)
フィルター： i′ バンド (中心波長 0.77µm)
露出時間： 180 秒
EXP ID
赤経 (h:m:s)
赤緯 (d:m:s)
黄経 (deg)
黄緯 (deg)
時刻（UT）
SUPE01354310
20:34:58.153
−21:44:01.87
305.77
−2.96
09:14:10
SUPE01354320
20:34:55.678
−21:43:14.43
305.77
−2.94
09:17:51
SUPE01354330
20:34:51.696
−21:43:32.49
305.75
−2.94
09:21:31
SUPE01354340
20:34:51.660
−21:44:31.07
305.75
−2.96
09:25:09
SUPE01354350
20:34:55.701
−21:44:49.38
305.76
−2.97
09:28:54
「対象天体」を見て分かるように、これらデータは太陽系研究のために取得されたものではない。しかし
以下のような条件を満たしており、小惑星検出にも適している。
• 黄道面に近い（黄緯 −3◦ ）
• 黄経が衝に近い（黄経差 ∼3◦ ）
• 銀河面から離れている（銀緯 −32◦ ）
• 適切なフィルター（i′ バンド）
• 観測時間（約 30 分間）、撮像回数（5 回）、撮像間隔（3–4 分）が適切
2.2
フラット
この観測では同日に dome フラットと twilight フラットの両方が取得されているが、今回は後者を使用
する。データの詳細は以下の通り。
2
2.3
EXP ID
時刻 (UT)
露出時間 (秒)
背景レベル
SUPE01354730
15:09:59
47.8
20972
SUPE01354740
15:11:32
37.1
20960
SUPE01354750
15:12:54
30.5
21405
SUPE01354760
15:14:16
24.4
21329
SUPE01354770
15:15:31
20.2
21575
SUPE01354780
15:16:38
17.2
22093
SUPE01354790
15:17:48
14.3
22320
SUPE01354800
15:18:46
12.2
22508
SUPE01354810
15:19:42
10.4
22535
SUPE01354820
15:20:39
8.9
22903
測光較正用データ
測光原点を決定するために背景の SDSS1 天体を参照する。しかし、ターゲット天体は SDSS 領域ではな
いため、その前後に同じフィルターで取得された SDSS 領域のデータを使用する。露出時間は全て 180
秒である。
EXP ID
天体
赤経 (h:m:s)
赤緯 (d:m:s)
時刻（UT)
Airmass
SUPE01354290
ABELL2261
17:22:29.679
+32:07:17.67
09:02:53
1.0989
SUPE01354370
ABELL2552
23:11:30.221
+03:35:19.29
09:37:21
1.7790
これらのデータは /mfst01b/teraity/SubaruSchool/data/ 下に置かれており、ds9 などの FITS ビュー
ワーで画像を表示させることができる。
例：
% ds9 -mosaic wcs /mfst01b/teraity/SubaruSchool/data/SUPA0135431?.fits
なお、各人の作業ディレクトリにこれらの生データをコピーする必要はない。代わりに、reduction と
いう名前のディレクトリを作成し、その下にデータのリストをコピーしておく。
% cd [directory for your work]
% mkdir reduction
% cd reduction
% cp /mfst01b/teraity/SubaruSchool/reduction/object raw.lis .
% cp /mfst01b/teraity/SubaruSchool/reduction/skyflat raw.lis .
% cp /mfst01b/teraity/SubaruSchool/reduction/photcal raw.lis .
1 Sloan
Digital Sky Survey <http://www.sdss.org/>
3
データ解析
3
3.1
一次処理
データ解析実習で行った手順に従って、Suprime-Cam 解析用ソフト「SDFRED2」を用いて画像データ
の整約処理を実行する。ただし本コースでは、画像の重ね合わせは行わず、
1. bias 差し引き（overscansub.csh）
2. 感度補正（ffield.csh）
3. 歪み補正（distcorr.csh）
4. sky 差し引き（skysb.csh）
の４段階のみでよい。
すでにデータ解析実習で一通りの手順を習得済みのはずなので、ここでは簡略化のために、一連の作業
を一括で行う処理スクリプトを使用する。以下の作業は reduction ディレクトリ下で行うこと。
［ twilight フラットの作成 ］
% python /mfst01b/teraity/SubaruSchool/script/skyflat.pyc skyflat raw.lis
⇒ sky mflat [DETECTOR].fits および sky mflat.lis が作成される
ファイル整理のため、以下のようなディレクトリを作って移動させる。その際、画像リストの内容も更
新しておく。
% mkdir skyflat
% mv sky mflat*.fits skyflat
% sed -i ’s/sky/skyflatY
=/sky/’ sky mflat.lis
［ ターゲット領域データの処理 ］
% python /mfst01b/teraity/SubaruSchool/script/sdfred.pyc object raw.lis
sky mflat.lis
⇒ fTo RH*****.fits, sgfTo RH*****.fits, object fTo.lis, object sgfTo.lis
% mkdir object fTo
% mv fTo RH*.fits object fTo
% sed -i ’s/fTo RH/object fToY
=/fTo RH/’ object fTo.lis
% mkdir object sgfTo
% mv sgfTo RH*.fits object sgfTo
% sed -i ’s/sgfTo RH/object sgfToY
=/sgfTo RH/’ object sgfTo.lis
4
［ 測光較正用データの処理 ］
% python /mfst01b/teraity/SubaruSchool/script/sdfred.pyc photcal raw.lis
sky mflat.lis
⇒ fTo RH****.fits, sgfTo RH****.fits, photcal fTo.lis, photcal sgfTo.lis
% mkdir photcal fTo
% mv fTo RH*.fits photcal fTo
% sed -i ’s/fTo RH/photcal fToY
=/fTo RH/’ photcal fTo.lis
% mkdir photcal sgfTo
% mv sgfTo RH*.fits photcal sgfTo
=/sgfTo RH/’ photcal sgfTo.lis
% sed -i ’s/sgfTo RH/photcal sgfToY
3.2
位置較正
一次処理後のターゲット領域画像（object sgfTo/sgfTo RH*****.fits）を ds9 で表示させ、既存の天体
カタログ（以下では Guide Star Catalog 2.22 を使用）を重ねてみよう。
% ds9 &
“File” → “Open” → Selection → “OK”
“zoom” → “fit”
“scale” → “zscale”
“Edit” → “Preferences” → “catalogs” → “Server”を “adac”に変更
“Analysis” → “Catalogs” → “Optical” → “GSC2.2”
プロットされたカタログ天体の位置が画像上の天体とずれているのが分かると思う。これは、FITS 画
像の画素座標（x, y ）を天球座標（赤道座標など） に関連付ける標準的な機構である World Coordinate
System (WCS) が正しく設定されていないためである。画素座標と天球座標の対応関係を較正し、FITS
ヘッダの情報を更新することにより、例えば (x, y) ↔ (RA, Dec) 間の相互変換が正確にできるように
なる。特に移動天体研究のためのデータ解析では、この作業は必須といってよい。
ここでは位置較正ソフトウエアの WCSTools（SAO）3 を使用して上記の作業を実行する。作業手順は
以下の通りである。
(1) 画像上の（比較的明るい）天体の位置を測定して座標リストを作成
(2) 位置標準星カタログから同じ天域の天体リストを取得
(3) ２つのリストを比較して同じ天体を対応付ける
(4) 画素座標と天球座標の変換係数を決定し、FITS ヘッダーに追記
手順 (1) には、光源カタログ生成ソフトウエアの SExtractor
4
を使用する。SExtractor の実行には、
各種設定を記述する「configuration ファイル（.sex）」と、出力カタログの内容を指定する「catalog
parameter ファイル（.param）」の２つが必要である。
2 http://gsss.stsci.edu/Catalogs/GSC/GSC2/GSC2.htm
3 http://tdc-www.harvard.edu/wcstools/
4 Source-Extractor
(Bertin, E. & Arnouts, S. 1996) <http://www.astromatic.net/software/sextractor>
5
configuration ファイル（ここでは wcscorr.sex という名前とする）は
% sex -d > wcscorr.sex
で作成することができる。ただし、これはデフォルトの設定になっているため、データや目的に応じて
適切な内容に書き換える必要がある。今回は以下の項目を次のように変更する（vi や emacs などのエ
ディタを使用するとよい）。
CATALOG NAME
wcscorr.cat # name of the output catalog
CATALOG TYPE
ASCII
# NONE,ASCII,ASCII HEAD, ASCII SKYCAT,
# ASCII VOTABLE, FITS 1.0 or FITS LDAC
PARAMETERS NAME
..
.
wcscorr.param
# name of the file containing catalog contents
DETECT MINAREA
30
# minimum number of pixels above threshold
DETECT THRESH
20
# <sigmas> or <threshold>,<ZP> in mag.arcsec-2
ANALYSIS THRESH
20
# <sigmas> or <threshold>,<ZP> in mag.arcsec-2
FILTER
N
# apply filter for detection (Y or N)?
40000.0
# level (in ADUs) at which arises saturation
..
.
SATUR LEVEL
特に重要なパラメータは DETECT THRESH と DETECT MINAREA の２つ。前者は検出閾値を示すパラメータ
で、バックグラウンドの RMS（root mean square）値の何倍以上のカウント値を持つ画素を天体からの
光と見なすか、を指定する。後者は、その閾値以上のカウント値を持つ画素が何個連なっていればそれら
を天体と見なすか、を指定するパラメータである。位置較正を行うためには、既存のカタログに載ってい
るような比較的明るい恒星（点光源）を抽出する必要があるので、DETECT THRESH と DETECT MINAREA
は大きめの値を設定するとよい。上記では参考値としてそれぞれ 20 と 30 が代入されているが、これで
うまく行かない場合は適宜調整してほしい。
一方、catalog parameter ファイル（wcscorr.sex の PARAMETERS NAME で指定。ここでは wcscorr.param
という名前とする）については、以下のような内容のファイルを（echo コマンドもしくはエディタで）
作成する。
XWIN IMAGE
YWIN IMAGE
MAG AUTO
FLAGS
これら 2 つのファイルを用意したら、SExtractor を実行する。
% sex object sgfTo/sgfTo RH120724object048 [DETECTOR].fits
-c wcscorr.sex
6
実行が成功すると、CATALOG NAME で指定した光源カタログ（ここでは wcscorr.cat）が作成される。
wcscorr.param で指定した通り、1、2 行目は光源重心の XY 座標、3 行目は自動開口測光で測定された
等級（ここでは相対値として使用）、4 行目は検出光源の flag が羅列されているはずである（各パラメー
タの詳細は SExtractor のマニュアル 5 を参照すること）。
検出された光源には位置較正には適さないものも含まれている。1 以上の flag が立っている光源は、何
かしら問題がある（複数の天体が混ざっている、飽和している、画像端にあって途切れている、など）
可能性が高い。また、等級の値が 99.0000 になっているものは正しく測光でされなかった光源を意味す
る。したがって、そのどちらかに当てはまる光源は除外する。加えて、WCSTools の位置較正コマンド
（imwcs）は最大でも 200 個の天体しか使わないので、リストを明るい光源順にソートし、上位 200 個の
みを選ぶ。なお、その後の作業では flag 値は不要なのでカットしておく。
% awk ’{if($4 == 0 && $3 < 99) printf "%8.3f %8.3f %8.4fY
=n",
$1, $2, $3}’ wcscorr.cat | sort -n -k 3,3 | head -200
> wcscorr.txt
リストがどのような光源を含んでいるのか、実際に確認してみよう。ds9 にプロットするために、下記
のコマンドで region ファイル 6 を作成する。
% awk ’{printf "CIRCLE %8.3f %8.3f 10Y
=n", $1, $2}’ wcscorr.txt
> wcscorr.reg
作成した wcscorr.reg を、先ほど起動した ds9 上で画像にプロットする。概ね適切そうな星が選ばれて
いれば問題ない。
“Frame” → “New Frame”
“File” → “Open” → Selection → “OK”
“zoom” → “fit”
“Region” → “Load Regions” → Selection (wcscorr.reg) → “OK”
これで手順 (1) が完了し、次に手順 (2)∼(4) に進む。下記のようにして WCSTools の位置較正コマンド
imwcs を実行する。WCSTools の最新版（ver. 3.9.27 ）が
/mfst01b/teraity/SubaruSchool/wcstools-3.9.2
に置いてあるので、これを使用する。
% /mfst01b/teraity/SubaruSchool/wcstools-3.9.2/bin/imwcs -d wcscorr.txt
-n 8 -w -h 200 -c gsc2 -o wsgfTo RH120724object048 [DETECTOR].fits
-q irst sgfTo RH120724object048 [DETECTOR].fits
各引数の詳細 8 は以下の通り。
5 https://www.astromatic.net/pubsvn/software/sextractor/trunk/doc/sextractor.pdf
6 http://ds9.si.edu/doc/ref/region.html
7 http://tdc-www.harvard.edu/software/wcstools/wcstools-3.9.2.tar.gz
8 http://tdc-www.harvard.edu/software/wcstools/imwcs/imwcs.com.html
7
-v
-d
-n
-w
-h
-c
-o
-q
Print more information about process
Use this catalog
Fit number parameters to a plane tangent WCS
Write new file or header
Maximum number of reference stars to use (10-200)
Reference catalog
Write output to the file in the next argument
Add optional processing steps:
i: iterate
r: recenter
s: sigma clip
t: tolerance reduce
「-c gsc2」は、比較カタログに Guide Star Catalog 2.3 の使用を指定している。コマンド実行後、位置較
正が完了すると「-o」で指定した名前の画像ファイルが生成される。この画像を ds9 で表示し、“Analysis”
→ “Catalogs” で GSC2.2 などの天体カタログをプロットしてみよう。天体カタログと画像上の星がぴっ
たり合っていれば成功である。
ここまで来ると sgfTo RH*****.fits は不要なので、ディレクトリを作って移動させておく。
% mkdir object sgfTo
% mv sgfTo RH*.fits object sgfTo
% sed -i ’s/sgfTo RH/object sgfToY
=/fTo RH/’ object sgfTo.lis
また、wsgfTo RH*****.fits は名前がやや長いので、簡潔な名称に変更しておこう。
% python /mfst01b/teraity/SubaruSchool/script/rename.pyc object
これで wsgfTo RH*****.fits は object01 0.fits などに改名される（“01 0” は 1 番目に撮られたフレー
ムの chip 0 を表す）。これらも同じくディレクトリを作って移動させておく。
% mkdir object
% mv object*.fits object
8
4
移動天体検出
移動天体を見つけるには、同じ領域を時間をおいて撮像した画像を比較し、位置が変化していく天体を
探せばよい。これには様々な手法があり、対象天体、観測方法、データ量、解析環境などによって最適
なものを選択する。最も簡単な方法は、画像を連続的に表示（ブリンク）させ、目で見て移動天体を探
すことだが、データが大規模化するとたちまち限界を迎えてしまう。この模擬研究で実際に扱うデータ
量は大変小さいものの、Hyper Suprime-Cam などによって取得される広域サーベイデータへの拡張を
念頭に置き、効率的かつ均質なサーベイが可能な、移動天体の自動抽出を核とする方法で小惑星の検出
を行う。主な手順は以下の通り。
(1) 全てのフレームに対して光源検出ソフトウエアを動作させ、光源リストを作成
(2) フレーム間で光源の位置を比較し、一致したものは背景天体（恒星・銀河など）と認定してリス
トから除外
(3) フレーム間で光源の組み合わせを作り、小惑星の移動速度に合致したものを候補天体として選定
(4) 候補天体を画像上で確認し、移動天体かどうかを判定
今回使用しているデータは、同じ天域を 5 回（以下、5 visits と表記）連続して撮像しているため、メイ
ンベルト小惑星の抽出は比較的容易に行うことができる。visit 数が少なかったり、撮像間隔が必要以上
に長くなるほど、抽出作業はより困難になる。ただし、visit 数を多くすると観測できる視野数が少なく
なり、検出総数が減ってしまう。また、研究対象が外縁天体のように移動速度の遅い（∼3 arcsec/hr）天
体であれば、撮像間隔を長く（∼1 hr 以上）しなければならない。このような観測を実施する際は、様々
な条件を勘案し、限られた観測時間の中で研究課題の要求を満たすデータを質・量ともに最大化できる
よう、十分に検討しなければならない。
9
4.1
光源リスト作成
3.2 章で行ったように、SExtractor を使って光源リストを作成する。ただし今回は、なるべく暗い光源
まで検出することに重点を置く。
% sex -d > wcscorr.sex
で configuration ファイルを生成し、以下の項目を修正する。
CATALOG NAME
detect.cat
# name of the output catalog
CATALOG TYPE
ASCII
# NONE,ASCII,ASCII HEAD, ASCII SKYCAT,
# ASCII VOTABLE, FITS 1.0 or FITS LDAC
ARAMETERS NAME
detect.param
# name of the file containing catalog contents
DETECT MINAREA
5
# minimum number of pixels above threshold
DETECT THRESH
2.5
# <sigmas> or <threshold>,<ZP> in mag.arcsec-2
ANALYSIS THRESH
2.5
# <sigmas> or <threshold>,<ZP> in mag.arcsec-2
FILTER
N
# apply filter for detection (Y or N)?
PHOT APERTURES
25
# MAG APER aperture diameter(s) in pixels
SATUR LEVEL
40000.0
# level (in ADUs) at which arises saturation
BACK FILTERSIZE
4
# Background filter: <size> or <width>, <height>
DETECT THRESH を小さくするほど暗い天体が検出されるが、背景ノイズも多く拾ってしまい、移
動天体検出を難しくしてしまう。概ね 2.0–3.0 の値が適当だと思われるが、必要があれば適宜調整して
ほしい。
次に「detect.param」というファイル名で catalog parameter ファイルを作成する。内容は以下の通り。
NUMBER
XWIN IMAGE
YWIN IMAGE
ALPHAWIN SKY
DELTAWIN SKY
FWHM IMAGE
ELLIPTICITY
FLAGS
FLUX APER
FLUXERR APER
detect.sex と detect.param を作成したら SExtractor を実行。なお、ここからの作業はチップごとに行うの
で、どのチップにするかを決めてから作業を進めること。以下では chip 2（fio）のデータ（object/object0[1-
5] 2.fits）を使った場合のコマンドが記されている。
% sex object/object01 2.fits -c detect.sex -CATALOG NAME detect01 2.cat
% sex object/object02 2.fits -c detect.sex -CATALOG NAME detect02 2.cat
% sex object/object03 2.fits -c detect.sex -CATALOG NAME detect03 2.cat
% sex object/object04 2.fits -c detect.sex -CATALOG NAME detect04 2.cat
% sex object/object05 2.fits -c detect.sex -CATALOG NAME detect05 2.cat
10
実行後には、コマンドライン中の -CATALOG NAME で指定した名前のカタログ（detect0[1-5] 2.cat）が作
成されているはずである。それらから region ファイルを作って ds9 上で画像にプロットし、検出された
光源を確認しておこう。
% awk ’{printf "CIRCLE %11.7fd %+11.7fd 2.0Y
=" Y
=n", $4, $5}’ detect01 2.cat
> detect01 2.cat.reg
% awk ’{printf "CIRCLE %11.7fd %+11.7fd 2.0Y
=" Y
=n", $4, $5}’ detect02 2.cat
> detect02 2.cat.reg
% awk ’{printf "CIRCLE %11.7fd %+11.7fd 2.0Y
=" Y
=n", $4, $5}’ detect03 2.cat
> detect03 2.cat.reg
% awk ’{printf "CIRCLE %11.7fd %+11.7fd 2.0Y
=" Y
=n", $4, $5}’ detect04 2.cat
> detect04 2.cat.reg
% awk ’{printf "CIRCLE %11.7fd %+11.7fd 2.0Y
=" Y
=n", $4, $5}’ detect05 2.cat
> detect05 2.cat.reg
4.2
移動天体候補の抽出
作成された全フレームの光源カタログから、手順 (2) および (3) を経て移動天体候補を選び出す。本研究
の肝となる重要な工程であるが、今回は簡略化のために予め用意された抽出プログラムを使用する（実
際の研究ではこれを自分たちで製作する必要がある）。まず、使用するデータ（下記の例では chip 2）の
リストを用意する。
% ls -1 object/object* 2.fits > object chip2.lis
次に、以下のコマンドを実行する。
% python /mfst01b/teraity/SubaruSchool/script/modetect.pyc object chip2.lis
処理が完了すると、候補光源のリスト（mocand 2.dat）およびそれらの座標を示す region ファイル
（moc0[1-5] 2.reg）が作成される。mocand 2.dat には、抽出された移動天体候補ごとに、1 行目に候補
光源 ID（e.g. MOC001）、2 行目以降に各 visit 画像上の位置情報が
visit 番号
画像ファイル名
X 座標
Y 座標
赤経
赤緯
黄経
黄緯
というフォーマットで記述されており、末行には黄経・黄緯方向の移動速度とその決定誤差が arcsec/hr
で書かれている。誤検出の数が多過ぎる場合には、detect.sex のパラメータを調整する必要がある。そ
のような場合は講師に相談してもらいたい。
なお、このプログラムは黄道面上の衝の位置でのメインベルト小惑星検出に特化した仕様になっており、
黄経方向の速度範囲： −60 < vλ < −20 arcsec/hr
黄緯方向の速度範囲： −50 < vβ < +50 arcsec/hr
に適合する移動天体のみを抽出する。また、5 visits 中最低 3 フレームからの検出が可能である。
11
4.3
目視判定
光源カタログから抽出された移動天体候補には、ノイズや明るい天体周りの光などに起因する誤検出も
少なからず含まれている。これらを機械的な方法で正確に判別するのはなかなか難しく、直接目で見て
確認するのが確実である。
ds9 で全 visit の画像（object/object0[1-5] 2.fits）を表示し、region ファイル（moc0[1-5] 2.reg）をプロッ
トしてから、
“Frame” → “First frame”
“Frame” → “Match” → “WCS”
“Frame” → “Lock” → “WCS”
とすると、全画像の表示範囲が 1 visit 目の画像と同じ場所・スケールに固定されるので、効率良く確認
作業を進めることができる。あとは “zoom” → “+”/“−” や “Edit” → “Pan” などのツールを活用して、
候補光源を順番に見ていってほしい。
慣れないうちは、移動天体かどうか迷うことが多々あるかもしれない。特に暗い天体の場合は判定が難
しい。判断に困った時は講師や他のメンバーに相談するようにしよう。
移動天体と確認できた候補光源については、 mocand 2.dat 中の記述をコピーして「mvo chip2.txt」と
いう名前のリストファイルに貼りつけていく。このようにして移動天体リストが作成される。
12
測定
5
この模擬研究の命題は、小惑星の直径と個数の関係を調べることである。恒星などと比べると圧倒的に
地球に近いとはいえ、小惑星は大変小さな天体なので、Suprime-Cam の解像度では空間分解されず、画
像上で天体サイズを直接測ることはできない。したがって、等級と距離から直径を推定することになる。
等級は、画像上で天体が拡がる範囲内の画素値の合計を測定（測光）した後、それを物理的な光度へ変
換（測光較正）することで得られる。一方、距離は直接測ることはできないので、小惑星の軌道から算
出しなければならない。
5.1
天体中心
detect.param に XWIN IMAGE、YWIN IMAGE （ALPHAWIN SKY、DELTAWIN SKY も同様）と指定した場合、
SExtractor は検出した光源 S の中心位置を画素値 Ii (xi , yi ) の重心として以下のように算出する。
∑
∑
Ii yi
Ii xi
x̄ =
i∈S
∑
,
ȳ =
Ii
i∈S
i∈S
∑
Ii
i∈S
恒星や楕円銀河のようなフラックスが中心集中している天体像に対しては、この方法で正確に中心を決
めることができる。ところが、移動天体は動いた分だけ像が伸びるため、ピーク位置が必ずしも像中心
と一致しない。そのため、SExtractor の算出した中心座標が本来の中心位置からずれている場合がある。
まずは中心座標を画像上で確認しよう。先ほど作成した mvo chip2.txt から region ファイルを作成する。
% grep object01 mvo chip2.txt | awk ’{printf "POINT %11.7fd %+11.7fd # POINT=CROSS
=n", $5, $6}’ > mvo01 2.reg
Y
······
ds9 でこれらの region ファイルを画像にプロットし、天体中心を捉えて
いるかを確認する。もし明らかに中心からずれている場合は、ds9 上で
“edit” → “region” を押して天体中心に合わせて円を描く（右図）。
円にカーソルを合わせてダブルクリックすると、円の中心座標および半
径を表示させることができる。
中心座標を決定したら、「meas mvo chip2.txt」という名前のファイルを作成し、
MVO–v[visit#]–c[chip#]–[obj#]
X
Y
（ 例： MVO–v01–c2-01
という書式で位置を追記していってほしい。
13
359.7
1228.4 ）
5.2
シーイング測定
恒星は本来、点光源であるはずだが、地球大気や望遠鏡・観測装置の光学系を通してやってきた光が最
終的に検出器によって受信される過程で、ある程度広がった像を持つことになる。これを点広がり関数
（PSF： point spread function）と呼び、その広がり度合いはシーイングとして評価される。
点光源の PSF は、一般的にガウス分布
( 2
)
1
x + y2
I(x, y) =
exp
−
2πσ 2
2σ 2
1
で表わされる。シーイングはその半値幅（FWHM：
Full-Width at Half-Maximum）で定義され、
√
FWHM = 2 2 ln 2σ ≈ 2.35σ
Ellipticity
0.8
0.6
0.4
0.2
02
と関係付けられる。
3
4
5
6
7
FWHM (pix)
測光の作業に入る前に、画像のシーイングサイズを把握しておく必要がある。いくつか方法があるが、
ここでは IRAF の psfmeasure というタスクを使用する。
シーイングサイズを測定するには、点光源（主に恒星）を選び出す必要がある。まずは画像上の光源の
様子を大まかに理解しておこう。先ほど作成した光源カタログ（e.g. detect01 2.cat）を利用して、検出
光源の PSF 分布を俯瞰する。
/mfst01b/teraity/SubaruSchool/script/ に置いてある fwhm-ellip.plt を作業ディレクトリにコピーし、
% sed "s/detect.cat/detect01 2.cat/g" fwhm-ellip.plt | gnuplot -persist
を実行すると、右上図のような画像が表示される。横軸は検出光源の FWHM、縦軸は楕円率（E ：
Ellipticity）9 を表す。図中の青い点線で囲ったところが点光源が集中していると思われる範囲である。
この結果から、FWHM < 4.0、E < 0.10 という範囲である程度点光源を絞り込めると考えられる。さ
らにシーイング測定には明るい天体が適しているため、以下のようにして候補光源リストを作成する。
% awk ’{if($6 < 4.0 && $7 < 0.10 && $8 == 0 && $9 > 100000) printf "%8.3f %8.3fY
=n",
$2, $3}’ detect01 2.cat > psc01 2.cat
次に IRAF を起動（初期設定および起動方法についての説明は省略）し、noao > obsutil というパッ
ケージをロードすると、psfmeasure タスクを使用することができるようになる。
obsutil> psfmeasure object/object01 2.fits coords="markall" imagecur="psc01 2.cat"
size="FWHM" radius=5 scale=1 saturation=35000 ignore sat+ displayというコマンドを実行すると、irafterm というウインドウが立ち上がり、X 軸（Column）および Y 軸
（Line）の FWHM 分布および Ellipticity 分布が表示される。それぞれ大半の点がある値付近に集中し
ているが、一部そこから外れている点もある。それらは点光源ではない可能性が高いので、その位置に
9 ellipticity
にはいくつか異なる定義がある。SExtractor については User’s manual の “Measurements” 章を参照。
14
カーソルを当てて “d” を押し、そのような点を除去する。それが完了したら、続いて “m” を押し、相
対等級 −FWHM ／ Ellipticity 分布を表示される。ここでも外れ値は “d” を押して除去しよう。さら
に “p” を押すと光源の radial profile を見ることができる。最後に “q” を押して処理を終了させると、
FWHM の平均値が表示される。この値をシーイングサイズと見なすことができる。
5.3
測光
開口測光（aperture photometry）：
The total number of counts attributable to the object is given by
Cobj =
n
∑
[ Cobj+sky (i) − Csky,e (i) ]
i=1
Cobj+sky (i)： number of counts in puxel i that are attributable
to the object and background (“sky”)
Csky,e (i)： estimated number of counts in the background of pixel i
( Reference: Newberry et al. 1991, PASP, 103, 122 )
移動天体の開口測光には、移動速度に応じて伸びた開口を使用するべきだが、今回は簡略化して点光源
用よりもやや大きな径の円形開口を採用する。また、同じく簡略化のため、SDFRED の skysb.csh 処
理によって全ての背景光が適切に差し引かれていると仮定し、Csky (i) = 0 とする。
開口径は主にシーイングサイズに合わせて決定するが、この解析では
点光源用の開口円半径：
15 pixels（3.0′′ ）
移動天体用の開口円半径：
22 pixels（4.4′′ ）
に統一する。
測光には IRAF の APPHOT パッケージにある phot タスクを使用する。パッケージを noao > digiphot
> apphot の順でロードした後、epar コマンドで各種パラメータを設定する。編集後、:wq とすると上
書き保存、 :q! とすると保存せずに終了する。
apphot> epar datapars
(readnoi=
(gain
=
10.)
GAIN)
CCD readout noise in electrons
CCD gain image header keyword
apphot> epar centerpars
(calgori=
none)
Centering algorithm
apphot> epar fitskypars
(salgori=
constant)
(skyvalu=
0.)
Sky fitting algorithm
User sky value
apphot> epar photpars
(apertur=
List of aperture radii in scale units
22.)
15
次に、測定する天体座標のリストファイルを作成する。
% grep v01 meas mvo chip2.txt | awk ’{printf "%6.1f %6.1fY
=n", $2, $3}’ > coord01 2.txt
······
phot タスクを実行し、続いて txdump タスクを用いて測定結果ファイルから指定した値を抜粋する。
apphot> phot object/object01 2.fits coords="coord01 2.txt" output="phot.txt" verboseverify- interactiveapphot> txdump "phot01 2.txt" "xinit,yinit,rapert,flux,merr,perr" yes headersxinit、yinit は指定した天体座標、rapert は開口半径、flux は測定されたフラックス値、merr は等
級誤差、perr はエラーメッセージを表す。測定結果を meas mvo chip2.txt の当該天体に
[rapart]
[flux]
[merr]
という書式で追記（コピペ）しておこう。
最後に、
="/4.4Y
="/g’ > amvo01 2.reg
% cat mvo01 2.reg | sed ’s/2.0Y
% sed "s/CIRCLE/POINT/g" amvo01 2.reg | sed "s/4.4Y
="/# POINT=CROSS/g" > cmvo01 2.reg
という region ファイルを作成し、ds9 で表示した画像にプロットする。perr が何らかのエラーメッセー
ジを出している、もしくは測光範囲に別の光源がある場合、正しく測光されていないことが疑われる。問
題ない場合は “ o ” を、そうでない場合は “ x ” を meas mvo chip2.txt の各行末に追記しておくように。
16
6
測光較正
測定されたフラックスから、その天体の等級は以下のように求められる。
mag obj = magzero − 2.5 log Cobj + 2.5 log texp
magzero は等級原点（magnitude zero point）、texp は露出時間（秒）である。この magzero を等級が既
知の天体を用いて算出し、天体の測光値を等級に換算できるようにするのが測光較正である。等級原点
は大気の吸収・散乱によって変化する。その減光量は主に光が通過した大気の厚み（airmass）に比例し、
∆mag = a(λ) × airmass で表される（ a ( λ ) は波長 λ での比例係数）。一般的に airmass は 1/ cos z
（z は天頂角）で近似される。係数 a(λ) を決定するには、同一のフィルターでターゲット領域と同じ高
度、あるいは２つの異なった高度での等級原点を求める必要がある。
この模擬研究で使用しているデータは i′ フィルターで取得されたものであり、この場合、SDSS 天体を
用いて較正を行うのが簡便である。しかし、ターゲット領域は SDSS フィールド外に位置しているため、
今回の場合は後者のケースに該当する。
本作業には、 photcal ディレクトリに置いてある制約処理済みデータ（ここでは photcal01 2.fits と
photcal02 2.fits）を使用する。最初に、これらのフレーム ID と airmass のリストを作っておく。IRAF
の hselect タスクを使用し、以下のようにしてヘッダー情報を抜き出す。
ecl> hselect photcal/photcal0? 2.fits FRAMEID,AIRMASS yes > magz photcal chip2.txt
次に、以下のような内容で SExtractor の configuretion ファイル（ source.sex ）と parameter ファイル
（ source.param ）を用意する。
source.sex：
CATALOG NAME
source.cat
# name of the output catalog
CATALOG TYPE
ASCII
# NONE,ASCII,ASCII HEAD, ASCII SKYCAT,
# ASCII VOTABLE, FITS 1.0 or FITS LDAC
ARAMETERS NAME
source.param
# name of the file containing catalog contents
DETECT MINAREA
5
# minimum number of pixels above threshold
DETECT THRESH
5.0
# <sigmas> or <threshold>,<ZP> in mag.arcsec-2
ANALYSIS THRESH
5.0
# <sigmas> or <threshold>,<ZP> in mag.arcsec-2
FILTER
N
# apply filter for detection (Y or N)?
BACK FILTERSIZE
4
# Background filter: <size> or <width>, <height>
VERBOSE TYPE
QUIET
# can be QUIET, NORMAL or FULL
source.param：
ALPHAWIN SKY
DELTAWIN SKY
FLAGS
17
そして、以下のスクリプトを実行する。
% python /mfst01b/teraity/SubaruSchool/script/calibstar.pyc photcal/photcal01 2.fits
このプログラムは
(1) SDSS のデータベースにアクセスして、領域内の SDSS 天体を検索
(2) 周囲に別の光源が無い、適当な明るさ（ i′ = 19.0–20.5 mag ）の恒星を抽出
(3) 該当した恒星の位置と i′ 等級 10 、および r′ − i′ カラーを記述したファイルを出力
作成された calibstar photcal02 2.txt には、等級較正に適していそうな恒星がリストアップされている。
これを fix calibstar photcal02 2.txt という名前で複製し、以降の確認作業で不適と判断される天体は削
除していく。
ds9 で画像を表示させ、同じく作成された calibstar photcal02 2.reg をプロットする。表示された円（半
径 3′′ ）内に明らかに別の光源がある場合、この天体は適さないと見なす。さらに、
“Scale”
>
“Scale parameters”
>
Limits
Low: 35000
High: 35000
として円内に飽和している画素の有無を確認し、あれば除外する。
確認が完了したら、
% awk ’printf "%6.1f %6.1fY
=n", $4, $5’ fix calibstar photcal02 2.txt > coord02 2.txt
で座標リストを作成しつつ、IRAF の測光関連パラメータを
apphot> epar datapars
(readnoi=
(gain
=
10.)
GAIN)
CCD readout noise in electrons
CCD gain image header keyword
apphot> epar centerpars
(calgori=
centroid)
(cbox
=
3.)
Centering algorithm
Centering box width in scale units
apphot> epar fitskypars
(salgori=
constant)
(skyvalu=
0.)
Sky fitting algorithm
User sky value
apphot> epar photpars
(apertur=
List of aperture radii in scale units
15.)
というように設定し、phot タスクを実行。
apphot> phot photcal02 2.fits coords="coord02 2.txt" output="phot02 2.txt" verboseverify- interactiveapphot> txdump "phot02 2.txt" "flux,merr" yes headers- > flux photcal02 2.txt
ここで注意しなければいけないのが、観測装置と標準フィルターシステムの波長感度特性の違い（color
term）である。すなわち、Suprime-Cam の i′ バンドと SDSS の i′ バンドでは、天体の赤さ／青さに
よって測定される等級が差が生じるため、これを補正する必要がある。
10 AB
等級システム（Fν = 3631 Jy を 0 等級と定義； Oke & Gunn 1983）
18
i′ バンドの場合、−0.4 < r′ − i′ < 2.5 の範囲では
i′SDSS − i′SCam ≈ 0.10 (r′SDSS − i′SDSS )
と近似される（Yagi et al. 2010）。
これを踏まえ、選ばれた恒星から等級原点を算出し、それらの平均値を計算する。
% paste fix calibstar photcal02 2.txt flux photcal02 2.txt | awk ’{if($7>-0.4 && $7<2.5
&& $8>0) printf "%6.3fY
=n", $6-0.1*$7+2.5*log($8)/log(10)}’ > magz photcal02 2.txt
% awk ’{m+=$1} END{print m/NR}’ magz photcal02 2.txt
このようにして photcal01 2.fits と photcal02 2.fits から等級原点を測定すると、係数 a ( i′ ) を得られ、
これを使ってターゲット領域の等級原点を求めることが可能になる。
19
軌道・天体直径の推定
7
7.1
軌道計算
検出された小惑星の移動速度から、円軌道（軌道離心率 = 0）を仮定して軌道要素を算出する。ここで
は観測領域が衝（太陽との黄経差 180◦ ）かつ黄道面（黄緯 0◦ ）上に位置すると近似する。
小惑星の軌道長半径を a、軌道傾斜角を i とすると、黄経・黄緯方向の移動速度（ vλ , vβ ）は、
vλ
=
√
√
µ
√
( cos i − a ) ,
a(a − 1)
√
vβ
=
µ
√
sin i
a(a − 1)
と表される。なお、µ = GM⊙ （G は重力定数、M⊙ は太陽質量）である。
v=
√
vλ2 + vβ2 として、これを解くと、
a =
√
]
1 [ 2
√
√
4 − 4 µ v v2 − 4 µ v2
µ
v
−
v
,
v
−
2
λ
λ
β
2 v2
tan i =
vλ +
|vβ |
√
µ/(a − 1)
【参考】モンテカルロ・シミュレーションによる軌道計算精度
［ 軌道長半径の推定誤差 ］
［ 軌道傾斜角の推定誤差 ］
範囲
平均値
標準偏差
範囲
平均値
標準偏差
2.1 – 2.5 au
+0.04 au
0.16 au
0 – 10◦
−0.2◦
2.6◦
2.5 – 2.9 au
+0.01 au
0.11 au
10 – 20◦
+0.4◦
5.1◦
2.9 – 3.3 au
+0.00 au
0.09 au
20 – 30◦
+4.4◦
6.2◦
20
7.2
小惑星直径
距離 r の天体の光度（単位時間あたりに放出される全エネルギー）を L 、フラックス（単位面積、単
位時間あたりに観測者が受けるエネルギー）を F とすると、
F =
L
4πr2
太陽 − 地球間の距離を r⊕ とすると、太陽のフラックスは
F⊙ =
L⊙
2
4πr⊕
日心距離 r、直径 D の小惑星に入射する単位時間あたりの太陽光エネルギーは
( )2
D
L⊙
L⊙ D2
Lin = π
=
,
2
2
4πr
16r2
全入射エネルギーに対する全反射エネルギーの割合（bond albedo）を A とすると、反射光の放射率は
Lout = ALin =
AL⊙ D2
16r2
距離 ∆ から観測される小惑星のフラックスは
F =
Lout
AL⊙ D2
pL⊙ D2
=
=
4π∆2
64πr2 ∆2
16πr2 ∆2
p は幾何学反射率（geometric albedo）で
p =
2
∫π
0
A
A
=
4
sin θ dθ
ここでは簡略化のため、位相角（太陽–天体–観測者のなす角）による減光効果は無視する。
太陽との等級差は
m − m⊙ = −2.5 log
2
p D 2 r⊕
F
= −2.5 log
F⊙
4 r 2 ∆2
小惑星の絶対等級（r = 1 au、∆ = 1 au、位相角 = 0◦ に規格化した等級）は
H
= m − 5 log r ∆
= m⊙ + 5 log 2 r⊕ − 2.5 log p − 5 log D
したがって、小惑星の直径は
log D = 0.2 m⊙ + log 2 r⊕ − 0.5 log p − 0.2 H
• i′ バンドでの太陽の AB 等級：
• 太陽 − 地球間の距離：
m⊙ = −27.04 mag
(Blanton & Roweis 2007)
r⊕ = 1 au ≈ 1.496 × 10 km
• 小惑星の平均的な幾何学反射率：
8
p ∼ 0.1
21
8
8.1
付録
Suprime-Cam の装置仕様
検出素子
Hamamatsu Photonics KK （2048×4096 画素）
CCD 数
10 (5×2，ギャップ幅 14–16′′ )
ピクセルスケール
0.20′′
視野
34′ × 27′
ゲイン
2.5–3.7 e− /ADU
読み出しノイズ
10 e−
飽和レベル
150,000 e−
［FITS ヘッダー情報］
http://smoka.nao.ac.jp/about/fits/Sample data/SuprimeCam/SupCam sample.header
［Suprime-Cam の CCD 配置］
6
7
2
1
0
chihiro
clarisse
fio
kiki
nausicaa
8
9
5
4
3
ponyo
san
satsuki
sheeta
sophie
22
9
赤道座標系と黄道座標系
赤道座標系： 地球の赤道面を基準とした座標系で、赤経（ RA ; α ）と赤緯（ Dec ; δ ）で表わされる。
黄道座標系： 地球の公転面（黄道面）を基準とした座標系で、黄経（ λ ）と黄緯（ β ）で表わされる。
両座標系は以下のように関係付けられる。
cos β cos λ
= cos δ sin α
cos β sin λ
= sin δ sin ϵ + cos δ sin α cos ϵ
sin β
= sin δ cos ϵ − cos δ sin α sin ϵ
ϵ は赤道面に対する黄道の傾斜角で、ϵ ≈ 23.44◦
これを解くと、
−1
(
− sin β sin ϵ + cos β cos λ cos ϵ
cos β cos λ
α
= tan
δ
= sin−1 (sin β cos ϵ + cos β sin λ sin ϵ)
)
または、
(
sin δ sin ϵ + cos δ sin α cos ϵ
cos δ cos α
λ
= tan−1
β
= sin−1 (sin δ cos ϵ − cos δ sin α sin ϵ)
)
90
Dec (degree)
60
30
Ecliptic
0
-30
Milky way
-60
-90
24
21
18
15
12
9
6
RA (hour)
赤道座標系における黄道と銀河面。
23
3
0
9.1
小惑星の軌道分布
参照： 軌道要素データベース MPCORB (Minor Planet Center)
http://www.minorplanetcenter.net/iau/MPCORB.html
既知小惑星の軌道長半径分布
50000
Number
40000
30000
20000
10000
01
2
3
4
5
Semi-major axis (au)
既知小惑星の軌道離心率・傾斜角分布
40000
60000
Number
Number
30000
20000
20000
10000
0
0.0
40000
0.1
0.2
0.3
0.4
0.5
0.6
0
Eccentricity
0
10
20
Inclination (deg)
24
30
40
10
移動速度計算
以下のような小惑星を観測した場合、どのような速度（秒角／時）で
天球上を移動するのか算出せよ。
◦ 太陽を中心とする、半径 3 天文単位の円形（離心率 0）軌道
◦ 地球軌道と同一平面上を公転
◦ 太陽・地球・小惑星が一直線上に並ぶ位置関係にある時
2 次元回転座標系の動径・角度方向の単位ベクトルを er , eθ とすると、
ėθ = −θ̇ er
ėr = θ̇ eθ ,
位置ベクトル r = r er の速度および加速度は
ṙ = ṙ er + rθ̇ eθ
r̈ = (r̈ − rθ˙2 ) er + (2ṙθ̇ + rθ̈) eθ
「円形（離心率 0）軌道」の場合、ṙ = 0、θ̈ = 0 であるから、
r̈ = −rθ˙2 er
µ = GM⊙ （G は重力定数、M⊙ は太陽質量）とすると、運動方程式は
−
µ
ėr
r2
θ̇
= −rθ˙2 er
√
=
µ/r3
したがって、
ṙ =
√
µ/r eθ
地球の位置ベクトル (r, θ, ϕ) を R = (1, 0, 0) とすると、
Ṙ = (0,
√
µ, 0)
「地球軌道と同一平面上を公転」し、「太陽・地球・小惑星が一直線上に並ぶ位置関係にある時」の小惑
星（軌道半径 a 天文単位とする）の位置と速度は
r = (a, 0, 0),
ṙ = (0,
25
√
µ/a, 0)
この時、地球−小惑星の距離は ∆ = a − 1 であるから、地球から見た小惑星の角速度は
ω
=
ṙ − Ṙ
(∆
=
0,
−
)
√ (
µ
1
1− √
,
a−1
a
)
0
√
√
すなわち、地球から観測される小惑星の移動速度は − µ (1 − 1/ a )/(a − 1)
√ √
√
計算しやすいように変形すると − µ / { a ( a + 1)} (a ̸= 1)
√
µ = |Ṙ| であるから、
360◦
√
≈ 148 arcsec/hr
µ =
1 yr
軌道が「半径 3 天文単位」の小惑星の移動速度は、
Motion along ecliptic latitude
(arcsec/hr)
148
−√ √
= −31.3 arcsec/hr
3 ( 3 + 1)
30 ◦
40
20 ◦
20
10 ◦
0◦
0
10 ◦
−20
20 ◦
−40
30 ◦
−60
2.1au
−50
2.4au
2.7au
−40
3.3au
3.0au
−30
Motion along ecliptic longitude (arcsec/hr)
上記の領域を観測した場合の円軌道小惑星の移動速度分布。
実線は軌道長半径ごと、破線は軌道傾斜角ごとの速度線。
26
© Copyright 2024 ExpyDoc
