OTO - Ohshima Tamashima Observatory- - guide_algorithm Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
ここでは、トランジット系外惑星の測光観測に特化した自動ガイドソフトを開発するためのアルゴリズムを考察する。
すなわち、その観測内容は、比較的短時間（数十秒から数分間）の露出時間の画像を、数時間（2時間から１晩）に渡り、延べ数１００フレームほど自動で撮像する場合である。この場合、目的天体の高度は、天頂から地平付近まで大きく変化する可能性がある。
!問題の所在
観測に使用する望遠鏡がシュミットカセグレン光学系の場合、望遠鏡の向きによって主鏡が光軸に対して傾くため、望遠鏡の鏡筒の機械的中心軸はうまく恒星を追尾しているつもりでも、画像上での目的天体は位置がずれていく。これは、精密な測光精度を要求される系外惑星の観測では、系統的な測光誤差を生じる原因となり、避けなければならない。ガイドエラーの許容量は、観測の期間を通して、せいぜい１０ピクセル程度の範囲に抑えたい。
この場合、観測用望遠鏡に同架したガイド望遠鏡によるオートガイドでは、二つの望遠鏡の光軸がずれるために機能しない。イド望遠鏡を使わず、測光観測用CCDカメラにガイド用チップを内蔵したSBIG社製のCCDカメラなら、うまくガイド星が確保できる場合はうまくいく可能性もあるが、系外惑星の測光観測の場合はわざとデフォーカスしているために、少し暗めのガイド星は見えなくなり、ガイドに使えないという問題点がある。
!観測画像そのものでガイディングする
そこで、測光観測している画像そのものを使いガイドすることができれば、いくら長時間の観測でも視野が移動することはなくなる。
!!前提
１フレームの露出時間内のガイディングは放棄し、望遠鏡の追尾精度に任せる。ミード製の望遠鏡では３０秒から２分間の露出でも結構流れることがあるが、あきらめて、測光アパーチャを大きく取ることで我慢する。
!!ガイド用ソフトに必要な機能
撮像完了直後の画像からガイド修正量を求め、望遠鏡に指示を出すために次の４つの機能が必要である。
!!!（１）星像検出と重心位置を計算する機能
画像から自動で星像を検出し、明るい方から３つの星の重心位置を求める。
下にRubyで実装したプログラムソース（ag05_1.rb）を載せる。
!!!（２）視野同定機能
テンプレートと問題の画像の２枚の視野を同定し、そのズレの量を検出する機能
!!!（３）ガイド用修正量を求める機能
（２）のズレ量を入力し、出力として、望遠鏡に与える最適な修正パラメータ量を求める機能
!!!（４）望遠鏡の操作
指定した入力値に応じて望遠鏡の向きを修正する機能
ここではTheSky6の機能を使い[[rubyから望遠鏡を制御する|operateLX200withTheSky6]]
!!!ガイド用テンプレート
観測開始ころの初期画像をファイル名として指定し、ガイド用テンプレート（参照用画像）として使用する。あるいは”template.fts”と言う名前にコピーしてもよい。
トランジット用画像では、星像は十分明るいので、極限等級まで検出する必要はなく、検出に使う閾値は画像全体の平均値レベル以上に設定してよいだろう。
!!雲の通過対策
雲が通過するなどして、星像検出がうまくできない場合は、へたに星を探さないで、望遠鏡の追尾に任せ、ガイドソフトとしてはじっとしている必要がある。
<<<
# ag05=1.rb
#FITSファイルから最も明るい星像を検出して重心を得る (C)O.ohshima2009.6.82009.9.2
#前提
#（１）画像はデフォーカスされている
#（２）ガイドの対象の星は、サチっていない範囲で画像中で一番明るい（ものから３つ）
#（３）ズレの１／２ずつ修正する（周期的誤差による振動防止）
# ここで使っているxy座標系は、通常のグラフでのもの（＝右がx+、上がy+）なので
#SBIGCCD座標系では、回る向きが逆になる
require 'fits'
require 'nimage'
class StarImage
def initialize(filename)
f=Fits.new
f.read(filename)
@im = f.data
@nx=f.data.shape[0]
@sorted_im=@im.sort
end
def get_g
#imaxで定めるカウント数より暗くて、一番明るい星像の重心を求める
#threshold_factorで星像検出の閾値を定める
# threshold = median + stddev*threshold_factor
#重心検出に成功した時の戻り値は、[x座標、y座標、総カウント数]
#総カウント数は、閾値より大きな値の総計
#失敗したら,falseが返る
imax=60000
threshold_factor=3
star_size=5
i=1
until @sorted_im[-1*i]<imax
i=i+1
end
d=@im.eq(@sorted_im[-1*i]).where[0]
j=(d/@nx).to_i
i=d-j*@nx
@thres=@im.median+threshold_factor*@im.stddev
#p @thres=image.mean
p self.is_star?(i,j,star_size)
end
def is_star?(x,y,size)
#星像であるか否か？
#引数（検証出発点のXY座標, 濃度検出閾値、星像サイズ閾値FWHM）
#星像サイズより大きければ星像
@lx, @ly=fwd2wall(x,y) #星像の壁まで突き当たり、左手法を開始する座標
# print "start left_hand_search at (",@lx+1,",",@ly+1,")","\n" #in ccd pix coordinates
@maxx=@minx=@lx
@maxy=@miny=@ly
@gx_sum=@gy_sum=@im_sum=0.0
left_hand_search(@lx,@ly+1,0,1) #開始地点から右回りに探査を開始
#p "\n", 'now start residual search'
# left_hand_search(@lx,@ly-1,0,-1) #開始地点の左ピクセルから残りの探査を開始
if ((@maxx-@minx)>size and (@maxy-@miny)>size) then
#print "minimax:", @minx,",",@maxx,",",@miny,",",@maxy, "\n" #連続領域範囲
(@minx..@maxx).each do |x| #連続領域のピクセル値の加算
(@miny..@maxy).each do |y|
if @im[x,y]<@thres then
@im[x,y]=0
else
@im[x,y]=@im[x,y]-@thres
end
@im_sum = @im_sum + @im[x,y]
@gx_sum = @gx_sum + x*@im[x,y]
@gy_sum = @gy_sum + y*@im[x,y]
end
end
gx = @gx_sum/@im_sum #星像の重心計算
gy = @gy_sum/@im_sum
return gx, gy, @im_sum #戻り値
else
false
end
end
def left_hand_search(x,y,dxp,dyp) #再帰を使った連続領域範囲の探査
# 引数：現在のx,y座標、dxp:直前のx増分、dxy:直前のy増分
# 事前に定義が必要なインスタンス変数
# @im[] 披検査画像（Narray）
# @thres 閾値
pixel_value=@im[x,y]
# print "("+x.to_s+","+y.to_s+")", pixel_value
if pixel_value>= @thres then
if @miny>y then @miny=y end #星像範囲設定
if @maxx<x then @maxx=x end
if @maxy<y then @maxy=y end
if @minx>x then @minx=x end
unless (@lx==x and @ly==y) then
l_value=@im[x+dyp,y-dxp]
fw_value=@im[x+dxp,y+dyp]
r_value=@im[x-dyp,y+dxp]
case
when l_value >= @thres
#p 'toL'
left_hand_search(x+dyp,y-dxp,dyp,-dxp)
when fw_value >= @thres
#p 'fwd'
left_hand_search(x+dxp,y+dyp,dxp,dyp)
when r_value >= @thres
#p 'toR'
left_hand_search(x-dyp,y+dxp,-dyp,dxp)
else
#p 'back'
left_hand_search(x-dxp,y-dyp,-dxp,-dyp)
end
end
else
#p 'less than threshold'
@lx=@lx-dxp
@ly=@ly-dyp
left_hand_search(x-dxp,y-dyp,-dxp,-dyp)
end
end
def fwd2wall(x,y) #内部の点から星像の壁までx座標を進める
pixel_value=@im[x+1,y]
if pixel_value>=@thres
fwd2wall(x+1,y)
else
return x, y #この戻り値から左手法を開始する
end
end
end
#t0=Time.now
im=StarImage.new('./data/HAT-P-2.00000030.FIT')
#im=StarImage.new('./data/Wasp14b.00000060.SAO_83398.FIT')
im.get_g
#p Time.now-t0
>>>
<<<
上のag05.rbの動作記録
（１）デバッグ用のプリントアウト機能を使い、途中経過を出力したもの
o2@am2:~/ruby/fits$ ruby ag05.rb
2879.1038131781
start left_hand_search at (428,113)
(427,113)2583.0"less than threshold"
(427,112)2901.0"toL"
(426,112)5009.0"toL"
(426,113)3400.0"toR"
(425,113)5745.0"toL"
(425,114)3278.0"toR"
(424,114)4780.0"fwd"
(423,114)6463.0"toL"
(423,115)3022.0"toR"
(422,115)3224.0"fwd"
(421,115)3137.0"fwd"
(420,115)3157.0"fwd"
(419,115)2908.0"toR"
(419,114)4661.0"toL"
(418,114)3247.0"toR"
(418,113)4483.0"toL"
(417,113)3006.0"toR"
(417,112)3392.0"toL"
(416,112)2892.0"back"
(417,112)3392.0"toL"
(417,111)3953.0"fwd"
(417,110)3933.0"fwd"
(417,109)3633.0"fwd"
(417,108)3064.0"toR"
(418,108)4400.0"toL"
(418,107)3436.0"toR"
(419,107)4166.0"fwd"
(420,107)5226.0"toL"
(420,106)3206.0"toR"
(421,106)3635.0"fwd"
(422,106)3467.0"fwd"
(423,106)3481.0"fwd"
(424,106)3024.0"toR"
(424,107)5201.0"toL"
(425,107)3905.0"toR"
(425,108)5841.0"toL"
(426,108)3618.0"fwd"
(427,108)2907.0"toR"
(427,109)3006.0"fwd"
(427,110)3182.0"fwd"
(427,111)3234.0minimax:416,427,106,115
[421.942153070307, 110.763682001345, 450463.348564148]
（２）最終結果のみを出力したもの
o2@am2:~/ruby/fits$ ruby ag05.rb
[421.942153070307, 110.763682001345, 450463.348564148]
>>>
----
[[技術的な覚書き]]へ戻る
[[OTO-Ohshima Tamashima Observatory-]]トップへ戻る
{{counter}}
すなわち、その観測内容は、比較的短時間（数十秒から数分間）の露出時間の画像を、数時間（2時間から１晩）に渡り、延べ数１００フレームほど自動で撮像する場合である。この場合、目的天体の高度は、天頂から地平付近まで大きく変化する可能性がある。
!問題の所在
観測に使用する望遠鏡がシュミットカセグレン光学系の場合、望遠鏡の向きによって主鏡が光軸に対して傾くため、望遠鏡の鏡筒の機械的中心軸はうまく恒星を追尾しているつもりでも、画像上での目的天体は位置がずれていく。これは、精密な測光精度を要求される系外惑星の観測では、系統的な測光誤差を生じる原因となり、避けなければならない。ガイドエラーの許容量は、観測の期間を通して、せいぜい１０ピクセル程度の範囲に抑えたい。
この場合、観測用望遠鏡に同架したガイド望遠鏡によるオートガイドでは、二つの望遠鏡の光軸がずれるために機能しない。イド望遠鏡を使わず、測光観測用CCDカメラにガイド用チップを内蔵したSBIG社製のCCDカメラなら、うまくガイド星が確保できる場合はうまくいく可能性もあるが、系外惑星の測光観測の場合はわざとデフォーカスしているために、少し暗めのガイド星は見えなくなり、ガイドに使えないという問題点がある。
!観測画像そのものでガイディングする
そこで、測光観測している画像そのものを使いガイドすることができれば、いくら長時間の観測でも視野が移動することはなくなる。
!!前提
１フレームの露出時間内のガイディングは放棄し、望遠鏡の追尾精度に任せる。ミード製の望遠鏡では３０秒から２分間の露出でも結構流れることがあるが、あきらめて、測光アパーチャを大きく取ることで我慢する。
!!ガイド用ソフトに必要な機能
撮像完了直後の画像からガイド修正量を求め、望遠鏡に指示を出すために次の４つの機能が必要である。
!!!（１）星像検出と重心位置を計算する機能
画像から自動で星像を検出し、明るい方から３つの星の重心位置を求める。
下にRubyで実装したプログラムソース（ag05_1.rb）を載せる。
!!!（２）視野同定機能
テンプレートと問題の画像の２枚の視野を同定し、そのズレの量を検出する機能
!!!（３）ガイド用修正量を求める機能
（２）のズレ量を入力し、出力として、望遠鏡に与える最適な修正パラメータ量を求める機能
!!!（４）望遠鏡の操作
指定した入力値に応じて望遠鏡の向きを修正する機能
ここではTheSky6の機能を使い[[rubyから望遠鏡を制御する|operateLX200withTheSky6]]
!!!ガイド用テンプレート
観測開始ころの初期画像をファイル名として指定し、ガイド用テンプレート（参照用画像）として使用する。あるいは”template.fts”と言う名前にコピーしてもよい。
トランジット用画像では、星像は十分明るいので、極限等級まで検出する必要はなく、検出に使う閾値は画像全体の平均値レベル以上に設定してよいだろう。
!!雲の通過対策
雲が通過するなどして、星像検出がうまくできない場合は、へたに星を探さないで、望遠鏡の追尾に任せ、ガイドソフトとしてはじっとしている必要がある。
<<<
# ag05=1.rb
#FITSファイルから最も明るい星像を検出して重心を得る (C)O.ohshima
#前提
#（１）画像はデフォーカスされている
#（２）ガイドの対象の星は、サチっていない範囲で画像中で一番明るい（ものから３つ）
#（３）ズレの１／２ずつ修正する（周期的誤差による振動防止）
# ここで使っているxy座標系は、通常のグラフでのもの（＝右がx+、上がy+）なので
#SBIGCCD座標系では、回る向きが逆になる
require 'fits'
require 'nimage'
class StarImage
def initialize(filename)
f=Fits.new
f.read(filename)
@im = f.data
@nx=f.data.shape[0]
@sorted_im=@im.sort
end
def get_g
#imaxで定めるカウント数より暗くて、一番明るい星像の重心を求める
#threshold_factorで星像検出の閾値を定める
# threshold = median + stddev*threshold_factor
#重心検出に成功した時の戻り値は、[x座標、y座標、総カウント数]
#総カウント数は、閾値より大きな値の総計
#失敗したら,falseが返る
imax=60000
threshold_factor=3
star_size=5
i=1
until @sorted_im[-1*i]<imax
i=i+1
end
d=@im.eq(@sorted_im[-1*i]).where[0]
j=(d/@nx).to_i
i=d-j*@nx
@thres=@im.median+threshold_factor*@im.stddev
#p @thres=image.mean
p self.is_star?(i,j,star_size)
end
def is_star?(x,y,size)
#星像であるか否か？
#引数（検証出発点のXY座標, 濃度検出閾値、星像サイズ閾値FWHM）
#星像サイズより大きければ星像
@lx, @ly=fwd2wall(x,y) #星像の壁まで突き当たり、左手法を開始する座標
# print "start left_hand_search at (",@lx+1,",",@ly+1,")","\n" #in ccd pix coordinates
@maxx=@minx=@lx
@maxy=@miny=@ly
@gx_sum=@gy_sum=@im_sum=0.0
left_hand_search(@lx,@ly+1,0,1) #開始地点から右回りに探査を開始
#p "\n", 'now start residual search'
# left_hand_search(@lx,@ly-1,0,-1) #開始地点の左ピクセルから残りの探査を開始
if ((@maxx-@minx)>size and (@maxy-@miny)>size) then
#print "minimax:", @minx,",",@maxx,",",@miny,",",@maxy, "\n" #連続領域範囲
(@minx..@maxx).each do |x| #連続領域のピクセル値の加算
(@miny..@maxy).each do |y|
if @im[x,y]<@thres then
@im[x,y]=0
else
@im[x,y]=@im[x,y]-@thres
end
@im_sum = @im_sum + @im[x,y]
@gx_sum = @gx_sum + x*@im[x,y]
@gy_sum = @gy_sum + y*@im[x,y]
end
end
gx = @gx_sum/@im_sum #星像の重心計算
gy = @gy_sum/@im_sum
return gx, gy, @im_sum #戻り値
else
false
end
end
def left_hand_search(x,y,dxp,dyp) #再帰を使った連続領域範囲の探査
# 引数：現在のx,y座標、dxp:直前のx増分、dxy:直前のy増分
# 事前に定義が必要なインスタンス変数
# @im[] 披検査画像（Narray）
# @thres 閾値
pixel_value=@im[x,y]
# print "("+x.to_s+","+y.to_s+")", pixel_value
if pixel_value>= @thres then
if @miny>y then @miny=y end #星像範囲設定
if @maxx<x then @maxx=x end
if @maxy<y then @maxy=y end
if @minx>x then @minx=x end
unless (@lx==x and @ly==y) then
l_value=@im[x+dyp,y-dxp]
fw_value=@im[x+dxp,y+dyp]
r_value=@im[x-dyp,y+dxp]
case
when l_value >= @thres
#p 'toL'
left_hand_search(x+dyp,y-dxp,dyp,-dxp)
when fw_value >= @thres
#p 'fwd'
left_hand_search(x+dxp,y+dyp,dxp,dyp)
when r_value >= @thres
#p 'toR'
left_hand_search(x-dyp,y+dxp,-dyp,dxp)
else
#p 'back'
left_hand_search(x-dxp,y-dyp,-dxp,-dyp)
end
end
else
#p 'less than threshold'
@lx=@lx-dxp
@ly=@ly-dyp
left_hand_search(x-dxp,y-dyp,-dxp,-dyp)
end
end
def fwd2wall(x,y) #内部の点から星像の壁までx座標を進める
pixel_value=@im[x+1,y]
if pixel_value>=@thres
fwd2wall(x+1,y)
else
return x, y #この戻り値から左手法を開始する
end
end
end
#t0=Time.now
im=StarImage.new('./data/HAT-P-2.00000030.FIT')
#im=StarImage.new('./data/Wasp14b.00000060.SAO_83398.FIT')
im.get_g
#p Time.now-t0
>>>
<<<
上のag05.rbの動作記録
（１）デバッグ用のプリントアウト機能を使い、途中経過を出力したもの
o2@am2:~/ruby/fits$ ruby ag05.rb
2879.1038131781
start left_hand_search at (428,113)
(427,113)2583.0"less than threshold"
(427,112)2901.0"toL"
(426,112)5009.0"toL"
(426,113)3400.0"toR"
(425,113)5745.0"toL"
(425,114)3278.0"toR"
(424,114)4780.0"fwd"
(423,114)6463.0"toL"
(423,115)3022.0"toR"
(422,115)3224.0"fwd"
(421,115)3137.0"fwd"
(420,115)3157.0"fwd"
(419,115)2908.0"toR"
(419,114)4661.0"toL"
(418,114)3247.0"toR"
(418,113)4483.0"toL"
(417,113)3006.0"toR"
(417,112)3392.0"toL"
(416,112)2892.0"back"
(417,112)3392.0"toL"
(417,111)3953.0"fwd"
(417,110)3933.0"fwd"
(417,109)3633.0"fwd"
(417,108)3064.0"toR"
(418,108)4400.0"toL"
(418,107)3436.0"toR"
(419,107)4166.0"fwd"
(420,107)5226.0"toL"
(420,106)3206.0"toR"
(421,106)3635.0"fwd"
(422,106)3467.0"fwd"
(423,106)3481.0"fwd"
(424,106)3024.0"toR"
(424,107)5201.0"toL"
(425,107)3905.0"toR"
(425,108)5841.0"toL"
(426,108)3618.0"fwd"
(427,108)2907.0"toR"
(427,109)3006.0"fwd"
(427,110)3182.0"fwd"
(427,111)3234.0minimax:416,427,106,115
[421.942153070307, 110.763682001345, 450463.348564148]
（２）最終結果のみを出力したもの
o2@am2:~/ruby/fits$ ruby ag05.rb
[421.942153070307, 110.763682001345, 450463.348564148]
>>>
----
[[技術的な覚書き]]へ戻る
[[OTO-Ohshima Tamashima Observatory-]]トップへ戻る
{{counter}}