※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

Difmapによる処理(5)CLEANによるdeconvolution


CLEANを実行する

CLEANによるdeconvolutionのアルゴリズムについて詳しく知りたい人は、こちらをご覧下さい。要するに、ダーティマップの中で輝度が最大の場所にCLEAN componentを置き、そのCLEAN componentから導かれるビジビリティを観測データから差し引いて残差ビジビリティにし、その残差ビジビリティを使って新たなダーティマップを描く…というプロセスを繰り返す反復法です。
さて、CLEANの効果を見やすくするために、ちょっと小細工をしてみましょう。残差マップをmapplotで表示するときに、輝度の値と画面上での明るさの対応は、デフォルトでは輝度の最小値・最大値を画面の最小輝度・最大輝度に対応させるように設定されています(この対応を transfer function: 伝達関数といいます)。この可変スケーリングは親切ですけど、CLEANによって残差が減っていくことを実感しづらいでしょう。そこで伝達関数を固定して、CLEANによって残差が減る様子を分かりやすくしてみます。それには、
mapfunc linear, -0.5, 3.0  第1引数は線型の伝達関数を指定, 第2引数は輝度の最小値, 第3引数は輝度の最大値を指定。
! Mapplot transfer-function = linear, Data range = -0.6 - 3 Jy.
と打ちます。これで輝度の伝達関数が線型に, 最小・最大輝度のレンジが-0.5 - 3.0 Jy/beam に設定されました。
この状態で、再度ダーティマップを表示しておきましょう。
mappl  
!
! Move the cursor into the plot window and press 'H' for help
PGPLOT画面上でxキーを押して、コマンドプロンプトに戻ります。それではいよいよcleanコマンドでCLEANします。
clean -1000, 0.001  CLEAN gainを0.001に設定して、1000回の反復でCLEANする。
! clean: niter=-1000 gain=0.001 cutoff=0
! Component: 050 - total flux cleaned = 0.166858 Jy
! Component: 100 - total flux cleaned = 0.325574 Jy
! Component: 150 - total flux cleaned = 0.476546 Jy
! Component: 200 - total flux cleaned = 0.620151 Jy
! Component: 250 - total flux cleaned = 0.756749 Jy
! Component: 300 - total flux cleaned = 0.886682 Jy
! Component: 350 - total flux cleaned = 1.01027 Jy
! Component: 400 - total flux cleaned = 1.12784 Jy
! Component: 450 - total flux cleaned = 1.23966 Jy
! Component: 500 - total flux cleaned = 1.34603 Jy
! Component: 550 - total flux cleaned = 1.44721 Jy
! Component: 600 - total flux cleaned = 1.54345 Jy
! Component: 650 - total flux cleaned = 1.635 Jy
! Component: 700 - total flux cleaned = 1.72208 Jy
! Component: 750 - total flux cleaned = 1.80491 Jy
! Component: 800 - total flux cleaned = 1.8837 Jy
! Component: 850 - total flux cleaned = 1.95864 Jy
! Component: 900 - total flux cleaned = 2.02993 Jy
! Component: 950 - total flux cleaned = 2.09774 Jy
! Component: 1000 - total flux cleaned = 2.16224 Jy
! Total flux subtracted in 1000 components = 2.16224 Jy
! Clean residual min=-0.261898 max=1.257721 Jy/beam
! Clean residual mean=0.000468 rms=0.100505 Jy/beam
! Combined flux in latest and established models = 2.16224 Jy
cleanコマンドの第1引数は反復回数です。マイナスを付けるのは、残差マップ中で最小輝度(負の値)の絶対値が最大輝度を上回ったらそこで反復を止める、というオプションです。つまり負のCLEAN成分を避けることができます。第2引数はCLEAN gainで、残差マップの最大値に対してCLEAN gain倍の値を、そこに置くCLEAN componentのフラックス密度とします。通常CLEAN gain は1より小さい値で、ここではとても慎重にCLEANを行うために0.001と小さな値にしています。CLEAN gainが大きすぎると、CLEANのしすぎ(CLEAN成分を差し引きすぎること)の原因になります。特に残差が大きいうち(CLEANがあまり進んでいない段階)ではCLEAN gainを控えめにした方がいいでしょう。強気に攻める場合でも0.1くらい、慎重な人は0.001くらいから始めてください。
さて、CLEAN実行時のメッセージを見てみましょう。total flux cleaned とあるのは、その段階までにCLEANした成分のフラックス密度の和です。理想的にはtotal flux cleanedが天体の全フラックス密度 (= 空間周波数ゼロλにおけるビジビリティ振幅 : radplでチェック!)になるまでCLEANを進めるべきでしょう。1000回の反復が終わった後で、残差マップの最小値が-0.26 Jy/beam, 最大値が1.26 Jy/beamですから、まだ正の側に偏っており、CLEANが十分ではないと判ります。

CLEANしたらチェック

このままCLEANを続けてもいいのですが、ここでCLEANの効果を見るためにmapplコマンドで残差を表示してみましょう。
mappl  
!
! Move the cursor into the plot window and press 'H' for help

中心の輝度がだいぶ小さくなり、残差が減っていることが分かるでしょう。ここでmキーを押すと、中心に緑色の十字が現れます。これがCLEAN componentです。「1000回も反復したのにCLEAN componentが1個しかないのはどうして!?」と思われるかも知れませんが、difmapでは同じ位置のCLEAN componentは1つにまとめてしまう仕様なのです。xキーでコマンドプロンプトに戻ってください。

edmod でCLEAN componentを見る

edmodコマンドで、CLEAN componentsをテキストエディターで扱うことができます。テキストエディターは、デフォルトでは vi です。difmapを起動する前に環境変数 EDITOR を設定しておけば、お好きなエディターを用いることができます。例えば emacs に慣れているのでしたら setenv EDITOR emacs とdifmap起動前に設定しておくとよいでしょう(起動した後に言うなよ!というツッコミは聞き流します)。それから、csh とか bash ではEDITOR="emacs"; export EDITORと設定してください。
ともあれ、edmodはCLEAN componentを閲覧したり編集したり保存したりすることができる便利なコマンドです。
! Flux (Jy) Radius (mas) Theta (deg) Major (mas) Axial ratio Phi (deg) T
! Freq (Hz) SpecIndex
2.0033 0.00000 0.00000
edmodで表示されるリストで、先頭が!の行はコメントです。したがって最後の行だけが意味があります。最初のカラム2.16224はCLEAN componentのフラックス密度、2・3番目のカラムで成分の重心の位置を表していて第2パラメーターはマップ中心からの距離 (mas), 第3パラメーターはその位置角(北を0$deg;として反時計回りに計る)です。第4 - 6パラメーターは、成分が楕円ガウシアンなど広がりを持っているときには意味がありますが、ここではCLEAN componentsをδ関数として扱っているので使いません。
CLEAN componentを確かめたら、エディターの終了コマンド (viなら :q) を打ってedmodを終了します。

さらにCLEANを進める

さらにCLEANを進めましょう。cleanと引数を省略すると、前回のCLEANで用いたパラメーターがそのまま引き継がれます。
clean
! clean: niter=-1000 gain=0.001 cutoff=0
! Component: 050 - total flux cleaned = 0.0613528 Jy
! Component: 100 - total flux cleaned = 0.119712 Jy
! Component: 150 - total flux cleaned = 0.175224 Jy
! Component: 200 - total flux cleaned = 0.228089 Jy
! Component: 250 - total flux cleaned = 0.278769 Jy
! Component: 300 - total flux cleaned = 0.327504 Jy
! Component: 350 - total flux cleaned = 0.374422 Jy
! Component: 400 - total flux cleaned = 0.419855 Jy
! Component: 450 - total flux cleaned = 0.463995 Jy
! Component: 500 - total flux cleaned = 0.506972 Jy
! Component: 550 - total flux cleaned = 0.548815 Jy
! Component: 600 - total flux cleaned = 0.589556 Jy
! Component: 650 - total flux cleaned = 0.629223 Jy
! Component: 700 - total flux cleaned = 0.667844 Jy
! Component: 750 - total flux cleaned = 0.705447 Jy
! Component: 800 - total flux cleaned = 0.742059 Jy
! Component: 850 - total flux cleaned = 0.777706 Jy
! Component: 900 - total flux cleaned = 0.812413 Jy
! Component: 950 - total flux cleaned = 0.846206 Jy
! Component: 1000 - total flux cleaned = 0.879108 Jy
! Total flux subtracted in 1000 components = 0.879108 Jy
! Clean residual min=-0.175274 max=0.649203 Jy/beam
! Clean residual mean=0.000386 rms=0.063478 Jy/beam
! Combined flux in latest and established models = 3.04134 Jy
CLEANしたらmapplotで確認。残差がさらに減っていることがわかりますね。mキーでCLEAN componentsを確認してみると、緑の十字がいくつか現れています。中心以外の場所にもCLEAN componentsが置かれたわけですね。確認したらxキーでコマンドプロンプトに戻ります。

だいぶ残差が小さく見づらくなってきましたから、そろそろ輝度伝達関数の自動スケーリングに戻しましょう。もう、CLEANの動作については分かったからいいですよね。
mapfunc linear, 0, 0  第2・第3引数を0にすると、自動スケーリング。
! Mapplot transfer-function = linear, Data range = data mim -> data max.

CLEAN boxを設定する

CLEANが進んで残差レベルが下がってきたら、雑音やサイドローブなど本来輝度が無いと思われるような場所にもCLEAN componentsを置かれる可能性が出てきます。それを避けるのが CLEAN boxです。CLEAN box は、囲った範囲の中にしかCLEAN componentsを置かないよう制限するためのものです。CLEAN boxは、mapplotの画面上で設定します。
mappl
!
! Move the cursor into the plot window and press 'H' for help
mapplotの画面上で、CLEAN boxの角に設定したい位置でマウス左ボタンをつんします。緑色の線がマウスポインターにまとわりつきますね。もう1つの角の位置で再度マウス左ボタンをつんすると、boxが設定できました。

CLEAN boxは複数設置しても構いません。重なり合っていても構いません。すでに存在するCLEAN boxを削除するには、消したいCLEAN boxの近くにマウスポインターを持って行ってDキーをタイプします。
CLEAN boxを設定する上で大切なことは、以下の点です。
  1. 「ここには輝度がないはず!」という場所を囲わないこと…そのためのCLEAN boxですから。
  2. 輝度が負の領域を囲わないこと…負のCLEAN componentが置かれてしまいます。
それでは、CLEAN boxを設定してCLEANを進めてください。
CLEANが煮詰まったら、Self Calibrationでさらにダイナミックレンジ向上を狙います。