Difmapによる処理(4) ダーティマップを描く
Dirty Mapを描く
いよいよマップを描きます。それにしても「ダーティ」ってどういうこと！？でしょう。もし空間周波数 (u, v) をすき間なくサンプリングできたなら完全なマップ（輝度分布I(l, m)）が描けるのですが、現実の観測では(u, v) coverageはすき間だらけで「完全」にはほど遠いです。空間周波数にすき間があると、合成ビームB(l, m)はサイドローブだらけです。このとき、得られるマップは真の輝度分布I(l, m)と合成ビームB(l, m)との畳み込み (convolution) になります。従って、得られるマップもサイドローブの影響で汚れます。これが「ダーティマップ」と言われるゆえんです。この辺の、真の輝度分布とダーティマップとの関係について詳しく知りたいと思ったら、
こちらのページをご覧下さい。
でも、ダーティマップからCLEANなどのアルゴリズムを使って真の輝度分布を推定することができます。これを deconvolution （デコンボリューション）といいます。convolutionを解くからdeconvolutionというのですね。CLEANによるdeconvolutionのアルゴリズムについて詳しく知りたい人は、
こちらをご覧下さい。
まあともあれ、Dirty Mapを描いて見てみましょう。それには前準備として以下の２つの手順が必要です。
- ビジビリティの重みづけを設定する (uvweight)
- マップの画素数とピクセルサイズを設定する (mapsize)
uvweightでビジビリティの重みづけを設定する
ビジビリティの逆フーリエ変換がダーティマップになるわけですが、手にしているビジビリティデータは離散的にサンプルされたビジビリティですので、離散逆フーリエ変換を行うわけです。その際に、各ビジビリティデータ点を全て同じ重み付けをするのが最適な結果を生む…とは限りません。ビジビリティは観測量ですから、それぞれのデータ点には誤差があります。この誤差が小さいデータには重みを大きく、誤差が小さいデータには重みを小さく与えるのが適切でしょう。統計的には、データ点の標準偏差σkに対して重み wk = 1/σk2を与えると、統計量としての標準偏差が最小にできます。ですが、重みが標準偏差の-2乗に比例するように設定すると、大小さまざまなアンテナが混在するアレイで観測したときには、小口径アンテナを含む基線のデータがほとんど使われないという悲しい事態におちいることがあります。ここではもう少しマイルドに、重みを標準偏差の-1乗に比例させる (wk ∝ σk-1) という方策をとってみます。
Difmapで重みづけを設定するには、uvweightコマンドを用います。
uvweigh 2, -1 2番目の引数が重みを標準偏差の-1乗に比例させることを指定。最初の引数は、正ならuniform weightを指定。
! Uniform weighting binwidth: 2 (pixels).
! Gridding weights will be scaled by errors raised to the power -1.
! Radial weighting is not currently selected.
uvweightの最初の引数の「2」は、(u, v)面のグリッドで2ピクセルの範囲で統計をとると指定しています（視野に関わるパラメーターなのですが、あまり深く追求しなくてもOK）。
mapsizeでマップの領域とピクセルサイズを指定
次にマップする領域の範囲とピクセルサイズを指定します。ピクセルサイズは、分解能（合成ビームサイズ）より十分に小さくとります。粗くてもビームFWHMの1/2以下 (Nyquist sampling)。できればFWHMの1/10 - 1/8程度が適切です。横（東西方向）と縦（南北方向）のピクセルサイズは違っていてもいいです。
マップの領域は、横（東西方向）と縦（南北方向）のピクセル数で指定します。ピクセル数は2nとなるような数でなくてはいけません。これは、フーリエ変換をFFT (高速フーリエ変換) で計算していることの都合です。縦と横のピクセル数は違っていてもいいです。
今回の観測では、基線長9000 kmのVLBAを使って波長2 cmで観測しているのですから、合成ビームのサイズは波長÷基線長で得られる分解能（ラジアン単位）をミリ秒角 (mas) 単位に直すことで、0.02 ÷ 9,000,000 = 2.2 ×
rad = 0.46 mas と見積もられます。そこで、ピクセルサイズを0.1 masに設定しましょう。
マップする領域を、256 × 256 ピクセル、すなわち22.6 × 25.6 masの領域に設定しましょう。観測天体のDA 193はコンパクトな電波銀河ですので、これくらいの狭い領域に十分収まります。ここで注意したいのは、マップするピクセル数の2倍をmapsizeコマンドに与えるということです。FFTの都合でmapsizeコマンドで与えた数の半分しか実際にマップに現れないのですが、ちょっと分かりにくいですよね。まあともかく、256ピクセルのマップ領域を得るためには、mapsizeコマンドで512と与える必要があります。
mapsize 512, 0.1 最初の引数はピクセル数の2倍, 2番目の引数がピクセルサイズ (mas)。
! Map grid = 512x512 pixels with 0.100x0.100 milli-arcsec cellsize.
もし縦と横でピクセル数やピクセルサイズを異なるように設定したいときには、横→縦の順で引数を与えます。例えば横はそのままで縦のピクセルサイズを0.2 mas, ピクセル数を512にしたいときには、mapsize 512, 0.1, 1024, 0.2とします。
mapplotでダーティマップを描く
mapplotを実行すると、ダーティマップが表示されます。
mappl
! Inverting map and beam
! Estimated beam: bmin=0.4196 mas, bmaj=1.028 mas, bpa=-3.999 degrees
! Estimated noise=1.21083 mJy/beam.
!
! Move the cursor into the plot window and press 'H' for help
メッセージからは、合成ビームの短軸が0.4196 mas, 長軸が1.028 mas, 長軸の位置角が-3.999°であることが読み取れます。また、雑音レベルがビジビリティの統計から1.21083 mJy/beamと見積もられることも分かります。これらの値は重要ですのでメモしておきましょう（ログファイルに自動的に記録されますけどね）。ここではダーティマップを表示させるだけで、CLEANすることもなく一旦コマンドプロンプトに戻りましょう。PGPLOT画面上で
xキーを押すと、コマンドプロンプトに戻ります。
mapplotで合成ビームを描く
mapplot実行時に beamというオプションを加えると、ダーティマップの代わりに合成ビームが表示されます。合成ビームのパターンを知っておくことで、どのような位置にサイドローブの影響が現れるかを把握しておきましょう。
mappl
mappl beam
!
! Move the cursor into the plot window and press 'H' for help
先ほどと同様に、PGPLOT画面上で
xキーを押して、コマンドプロンプトに戻ります。
ひきつづき、CLEANによるdeconvolutionを楽しみましょー！
最終更新：2007年10月04日 20:54
