イメージング

このチュートリアルの目的

このチュートリアルでは、CASA を初めて使う人向けに、ALMA データのイメージングの初歩を解説しています。この解説は、CASA guides "First Look at Imaging CASA 6.4" (本家のサイト:英語) の日本語訳(と多少の追加説明)となっています。

ALMA では、観測でデータが取得された後、パイプラインや ARC スタッフの作業によってキャリブレーション (大気や装置の影響を取り除く作業) が行われます。 そのキャリブレーションで使われたスクリプトは、ALMA アーカイヴ上で生データと一緒に配布されますので、ユーザーはキャリブレーション済みデータをすぐに入手できます。 ALMA の品質保証 (Quality Assurance 2 (QA2)) では、ビームサイズや到達感度の確認のためにイメージングも行われますが、サイエンスの目的に応じて、自分でイメージングをやり直したいと思うことも多いでしょう。 また、イメージングで用いられる CLEAN は、ALMA データに限らず、干渉計データの解析で頻繁に用いられる手法です。 このページと、他の一連のページ (セルフキャリブレーションラインのイメージングイメージの解析) で、イメージングの基本的な流れを示していきます。 なお、解析における基本的な概念、開口合成やデコンボリューションの解説、また、tclean 実行時に出てくるマスク、clean component、残差 (residual) などの意味については、例えば干渉計サマースクール2005 教科書Imaging Interferometric Data (英語)、その他、リンク集にある資料を参照してください。

CASA のインストール

まだ CASA をインストールしていなければ、こちらのページを参考にダウンロード、インストールしてください。このチュートリアルはバージョン 6.4 を対象として書かれていますので、そのバージョンをインストールして下さい。バージョンが多少違って動作する可能性は高いですが、多少見た目やオプションが違う可能性があります。

チュートリアルで用いるデータについて

例として用いるデータは、原始惑星系円盤を持つ若い星、うみへび座 TW 星 (TW Hydra) です。プロジェクトコード 2011.0.00340.S の "Searching for H2D+ in the disk of TW Hya v1.5" で得られたものです。2012 年 11 月にサブミリ波の連続光放射、H2D+ と N2H+ のラインをターゲットとして Band 7 で観測されました。 ここでイメージングするデータのバンド幅は 234.375 MHz で、610 kHz ごとに並んだ 384 個のチャンネルを含んでおり、12-m アンテナ 12 台を用いて Cycle 0 観測で得られたデータです。

サンプルデータのダウンロード

チュートリアル用のパッケージ (4.1 GB) は以下の場所から取得できます。
https://bulk.cv.nrao.edu/almadata/public/working/sis14_twhya.tgz

パッケージは tar コマンドで解凍できます。

# In bash
tar -xvzf sis14_twhya.tgz

このイメージングのチュートリアルでは、キャリブレーション後のデータを使います。上記を実行すると、"sis14_twhya_calibrated_flagged.ms.tar" というファイルができるので、これを更に展開し、データ "sis14_twhya_calibrated_flagged.ms" (measurement set, MS) を作業ディレクトリにコピーします。 そのディレクトリから CASA を起動するようにします (CASA ログは起動したディレクトリに保存されます)。

# In bash
cd working
tar -xvf sis14_twhya_calibrated_flagged.ms.tar
cd ..
mkdir MyTutorial
cd MyTutorial
cp -r ../working/sis14_twhya_calibrated_flagged.ms .

本ページに書かれているイメージングチュートリアル部分だけのデータ (600 MB) が欲しい場合は以下のように部分的に取得することもできます。

# In bash
mkdir MyTutorial
cd MyTutorial
wget -r -np -nH --cut-dirs=4 --reject "index.html*" https://bulk.cv.nrao.edu/almadata/public/working/sis14_twhya_calibrated_flagged.ms.tar

CASA を起動する

CASA を起動する前に、ls コマンドで sis14_twhya_calibrated_flagged.ms がカレントディレクトリにあるかどうかを確認しましょう。

# In bash
ls

CASA を起動します。

# In bash
casa

CASA を起動すると、CASA のロガー・ウィンドウが新たに立ち上がります。また、casa コマンドを打ったターミナルには CASA プロンプトが表示されます。なお、CASA は Python で動いていますので、Python の文法が使えます。バージョン 6 系は Python 3 になります。

CASA での作業を終えて CASA を終了するときは、exit() 又は quit() と入力します。

CASA の基本操作

CASA の操作に関するドキュメントとしては、CASA Guides の他の様々なページや cookbook があります (リンク集も参考にしてください)。ここでは、CASA を使い始めるにあたって知っておいた方が良いことを、少しだけ紹介します。 一部の Unix コマンドは、CASA でも利用できます。例えば、lspwd などです。

# In CASA
ls
pwd

次のように os.system を用いれば、すべての Unix/Linux コマンドを使うことができます。

# In CASA
os.system('ls')

また、CASA での解析操作を Python スクリプトとして記述しておき、それを一度に実行することもできます。スクリプトは、次のようにして実行します。

# In CASA
execfile('my_script.py',globals())

このチュートリアルは CASA ver. 6.4 を前提として書かれていますので、以下を実行して、CASA のバージョンを確認しておきます。

# In CASA
import casalith
version = casalith.version_string()
print ("You are using {}".format(version))
if (version < '6.4'):
    print ("YOUR VERSION OF CASA IS TOO OLD FOR THIS GUIDE.")
    print ("PLEASE UPDATE IT BEFORE PROCEEDING.")
else:
    print ("Your version of CASA is appropriate for this guide.")

タスクの使用とデータの確認

データ解析の手始めはデータのヘッダー情報やデータ構造を調べることです。listobs タスクを使用します。

CASA でタスクを実行する方法は 2 つあります。 パラメータを inp で予め設定しておき go で実行させる方法と、直接パラメータ込みのコマンドを実行する方法です。

# In CASA
inp listobs
vis='sis14_twhya_calibrated_flagged.ms'
go

# In CASA
listobs(vis='sis14_twhya_calibrated_flagged.ms')

タスクの使い方やパラメータがわからないときは、help で調べます。

# In CASA
help listobs

このヘルプは "q" で抜けることができます。 listobs の出力は以下のようになります。


================================================================================
                   MeasurementSet Name:  /lustre/naasc/sciops/qa2/dkim/casa_testing/MyTutorial/sis14_twhya_calibrated_flagged.ms      MS Version 2
================================================================================
   Observer: cqi     Project: uid://A002/X327408/X6f  
Observation: ALMA
Computing scan and subscan properties...
Data records: 80563       Total elapsed time = 5647.68 seconds
           Observed from   19-Nov-2012/07:36:57.0   to   19-Nov-2012/09:11:04.7 (UTC)

           ObservationID = 0         ArrayID = 0
          Date        Timerange (UTC)          Scan  FldId FieldName             nRows     SpwIds   Average Interval(s)    ScanIntent
          19-Nov-2012/07:36:57.0 - 07:39:13.1     4      0 J0522-364                 4200  [0]  [6.05] [CALIBRATE_BANDPASS#ON_SOURCE,CALIBRATE_PHASE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE]
                      07:44:45.2 - 07:47:01.2     7      2 Ceres                     3800  [0]  [6.05] [CALIBRATE_AMPLI#ON_SOURCE,CALIBRATE_PHASE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE]
                      07:52:42.0 - 07:53:47.6    10      3 J1037-295                 1900  [0]  [6.05] [CALIBRATE_PHASE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE]
                      07:56:23.5 - 08:02:11.3    12      5 TW Hya                    8514  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
                      08:04:36.3 - 08:05:41.9    14      3 J1037-295                 1900  [0]  [6.05] [CALIBRATE_PHASE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE]
                      08:08:09.6 - 08:13:57.3    16      5 TW Hya                   10360  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
                      08:16:20.6 - 08:17:26.2    18      3 J1037-295                 2100  [0]  [6.05] [CALIBRATE_PHASE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE]
                      08:19:53.9 - 08:25:41.7    20      5 TW Hya                   10321  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
                      08:28:17.1 - 08:29:22.6    22      3 J1037-295                 2100  [0]  [6.05] [CALIBRATE_PHASE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE]
                      08:32:00.5 - 08:37:48.2    24      5 TW Hya                   10324  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
                      08:40:11.9 - 08:41:17.4    26      3 J1037-295                 2100  [0]  [6.05] [CALIBRATE_PHASE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE]
                      08:43:45.6 - 08:49:33.4    28      5 TW Hya                    9462  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
                      08:51:57.1 - 08:53:02.6    30      3 J1037-295                 1900  [0]  [6.05] [CALIBRATE_PHASE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE]
                      08:58:12.0 - 09:00:28.1    33      6 3c279                     3402  [0]  [6.05] [CALIBRATE_BANDPASS#ON_SOURCE,CALIBRATE_PHASE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE]
                      09:01:35.7 - 09:02:41.2    34      3 J1037-295                 1900  [0]  [6.05] [CALIBRATE_PHASE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE]
                      09:05:15.6 - 09:07:31.6    36      5 TW Hya                    4180  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
                      09:09:59.1 - 09:11:04.7    38      3 J1037-295                 2100  [0]  [6.05] [CALIBRATE_PHASE#ON_SOURCE,CALIBRATE_WVR#ON_SOURCE]
           (nRows = Total number of rows per scan) 
Fields: 5
          ID   Code Name                RA               Decl           Epoch   SrcId      nRows
          0    none J0522-364           05:22:57.984648 -36.27.30.85128 J2000   0           4200
          2    none Ceres               06:10:15.950590 +23.22.06.90668 J2000   2           3800
          3    none J1037-295           10:37:16.079736 -29.34.02.81316 J2000   3          16000
          5    none TW Hya              11:01:51.796000 -34.42.17.36600 J2000   4          53161
          6    none 3c279               12:56:11.166576 -05.47.21.52464 J2000   5           3402
Spectral Windows:  (1 unique spectral windows and 1 unique polarization setups)
          SpwID  Name                           #Chans   Frame   Ch0(MHz)  ChanWid(kHz)  TotBW(kHz) CtrFreq(MHz) BBC Num  Corrs  
          0      ALMA_RB_07#BB_2#SW-01#FULL_RES    384   TOPO  372533.086       610.352    234375.0 372649.9688        2  XX  YY
Sources: 5 
          ID   Name                SpwId RestFreq(MHz)  SysVel(km/s)
          0    J0522-364           0     -              -            
          1    Ceres               0     -              -           
          2    J1037-295           0     -              -            
          3    TW Hya              0     -              -            
          4    3c279               0     -              -           
Antennas: 21:
          ID   Name  Station   Diam.    Long.         Lat.                Offset from array center (m)                ITRF Geocentric coordinates (m) 
                                                                             East         North     Elevation               x               y               z
          1    DA42  A050      12.0 m   -067.45.16.2  -22.53.29.3         43.0352     -744.9713       21.6702  2225079.880016 -5440041.377534 -2481724.598031  
          2    DA44  A068      12.0 m   -067.45.20.6  -22.53.25.7        -82.4232     -631.7828       23.5810  2224981.097784 -5440131.250387 -2481621.066374  
          3    DA45  A070      12.0 m   -067.45.11.9  -22.53.29.3        166.1833     -743.4934       19.8811  2225193.450167 -5439993.764157 -2481722.540534  
          4    DA46  A067      12.0 m   -067.45.12.7  -22.53.27.2        142.4097     -678.7318       20.1280  2225181.070532 -5440026.290790 -2481662.975103  
          5    DA48  A046      12.0 m   -067.45.17.0  -22.53.29.3         21.4267     -742.7987       21.6757  2225060.202580 -5440050.344436 -2481722.598651  
          6    DA49  A029      12.0 m   -067.45.18.2  -22.53.25.8        -12.9134     -636.4552       22.1350  2225044.239583 -5440102.022535 -2481624.808405  
          7    DA50  A045      12.0 m   -067.45.17.9  -22.53.30.1         -5.4183     -767.4398       22.6034  2225032.051652 -5440052.426015 -2481745.660003  
          9    DV02  A077      12.0 m   -067.45.10.1  -22.53.25.9        217.6299     -637.5333       15.8376  2225255.259272 -5440008.987869 -2481623.352052  
          11   DV05  A082      12.0 m   -067.45.08.3  -22.53.29.2        269.0433     -740.9521       15.7832  2225287.593766 -5439952.243679 -2481718.605314  
          12   DV06  A037      12.0 m   -067.45.17.5  -22.53.28.8          6.7403     -727.3003       21.2086  2225048.729287 -5440061.085777 -2481708.139136  
          14   DV08  A021      12.0 m   -067.45.17.2  -22.53.27.0         14.3196     -672.8108       21.3420  2225063.814715 -5440077.948261 -2481657.992572  
          15   DV10  A071      12.0 m   -067.45.19.9  -22.53.23.5        -60.7887     -563.2541       23.3799  2225011.141945 -5440147.560932 -2481557.855663  
          16   DV13  A072      12.0 m   -067.45.12.6  -22.53.24.0        147.1742     -580.5887       18.1825  2225199.254375 -5440058.161494 -2481571.803699  
          17   DV15  A074      12.0 m   -067.45.12.1  -22.53.32.0        161.8159     -828.6196       18.7688  2225176.483514 -5439963.820451 -2481800.529842  
          18   DV16  A069      12.0 m   -067.45.21.3  -22.53.30.2       -101.4797     -770.1047       23.2972  2224942.993176 -5440088.421459 -2481748.384855  
          19   DV17  A138      12.0 m   -067.45.17.1  -22.53.34.4         19.1461     -901.2603       26.0137  2225036.269025 -5439997.853009 -2481870.267607  
          20   DV18  A053      12.0 m   -067.45.17.3  -22.53.31.2         12.5939     -802.9941       21.5281  2225043.111690 -5440031.889497 -2481777.995870  
          21   DV19  A008      12.0 m   -067.45.15.4  -22.53.26.8         67.5592     -667.6872       20.9574  2225113.709955 -5440059.310545 -2481653.122797  
          22   DV20  A020      12.0 m   -067.45.17.8  -22.53.28.0         -2.9649     -703.4389       21.6629  2225043.419055 -5440073.737929 -2481686.333574  
          24   DV22  A011      12.0 m   -067.45.14.4  -22.53.28.4         95.9131     -716.5005       21.0898  2225132.810230 -5440031.115405 -2481698.143589  
          25   DV23  A007      12.0 m   -067.45.15.1  -22.53.27.3         74.0152     -681.2926       21.3231  2225117.809276 -5440052.280005 -2481665.79904

UV データをグラフィカルに確認するタスクは、plotms です。まず、データの UV カバレッジをプロットしてみます。

# In CASA
plotms(vis='sis14_twhya_calibrated_flagged.ms', xaxis='u', yaxis='v', avgchannel='10000', avgspw=False, avgtime='1e9', avgscan=False, coloraxis="field", showgui=True)

図 1 のように GUI が表示されます。UV カバレッジのプロットの色は、異なる観測ターゲットを示しています。

図 1: TW Hya データの UV カバレッジ

他にも様々なプロットを見ることができます。例えば、Amplitude (振幅) vs. UV distance を見るためには、"Axes" タブを選択し、"X Axis" を "UVDist" に、"Y Axis" を "Amp" に設定します。 "Page" タブで、"Iteration Axis" を "Field" に設定します。そして、Plot ボタンを押すと、プロットが表示されます。
plotms GUI の下部にある緑色の矢印ボタンを使って、様々なキャリブレーター (較正天体) とサイエンスターゲットに進むことができます。 ケレス (Celes) を含むいくつかの天体では、図 2 のように UV 距離の関数として振幅の変化を見ることができます。これらの天体は分解されており、下でおこなうイメージングではっきりと見ることができます。

図 2: ケレスの Amplitude 対 UV Distance のプロット

tclean を用いたイメージング

CASA 4.7 から、clean のイメージング機能は、tclean というタスクに改良されて置き換えられました。このガイドでは tclean のみを使用します。 TCLEAN と ALMA のガイド (英語) では、cleantclean の構文の違いが詳細に書かれています。 参考までに、イメージングに tclean ではなく clean を使用する古いバージョンの CASA Guide は、 Archived ALMA CASAguides に保管されています (英語のみ)。

このデータセットでは、位相較正天体はフィールド 3 で識別されます。この較正天体を "phase_cal" という画像ファイルにイメージ化することにします。 まず、以前このタスクを実行したことがある場合に備えて、この名前を使った古いバージョンの画像を削除しておきます。(画像化すると、拡張子が違うだけの同じルート名のファイルがいくつかできるので、".*" が必要です)。

# In CASA
os.system('rm -rf phase_cal.*')

それでは、tclean タスクを使ってイメージングをおこないます。tclean オプションの完全なリストはドキュメントを参照して下さい。
  • 最初に位相較正天体をイメージングし、イメージ名を "phase_cal" に設定します。MS 内の位相較正天体を特定するには、listobs の出力の "ScanIntent" 列を見ます。 位相較正天体は、"CALIBRATE_PHASE#ON_SOURCE" で示されています。この MS では、位相較正天体は J1037-295 であり、field = "3" で識別されます。
  • 多周波合成 (specmode = mfs) を使用して、単一の連続波イメージを作成します。多周波合成は、選択されたすべてのスペクトルチャンネルからのデータを 1 つの連続波イメージに合成します。 観測が広い周波数範囲をカバーする場合、周波数によって天体の振幅や構造が大きく変化する可能性があります。これは、バンド幅の割合 (Δν/ν_center) が 10 % より大きい場合にのみ考慮する必要があります。 このデータセットでは問題にはなりません。そのため、deconvolver = hogbom に設定し、nterms = 1 を使用し、tclean にデコンボリューションされた各成分がすべての周波数で単一の振幅を持つように設定します。
  • ポインティングは 1 つだけなので、gridder='standard' とします。複数の連続したポインティングがある場合、あるいは 7 [m] と 12 [m] のアレイデータを同時に撮像する場合 (それぞれ 1 ポインティングだけであっても)、gridder='mosaic' に設定する必要があります。
  • セルサイズは 0.1 [arcsec] とし、ビーム全体で 5 ピクセルとした。経験則として、楕円ビームの最小方向で 5 ピクセル程度が望ましいということが知られているためです。 小規模なデータセットの場合、クイックなイメージングで適切なセルサイズを決定し、tclean で計算されたビームを記録することができるかもしれません。 しかし、多くの ALMA プロジェクトでは、これにはかなりの時間がかかります。uv-coverage をもう一度見ることで、セルサイズを推定することができます。 x 軸を "uvwave" に変えると、セルサイズはおおよそ 206265/(波長で最も長いベースライン)/(ビーム全体のセル数) となります。 この MS の場合、これは 0.09 [arcsec] になり、0.1 [arcsec] に切り上げることにします。tclean は 0.58 [arcsec] × 0.51 [arcsec] の合成ビームサイズを出力しており、このセルサイズは適当であると言えるでしょう。
  • 画像サイズは 128 × 128 に設定しました (ただし CASA では 2 の累乗は特殊な数ではないので、別の大きさでも構いません)。 これは観測周波数における主ビームの大部分をカバーするのに十分ですが、非点源に対してはもっと広い視野を撮像したいかもしれません。画像サイズは任意ですが、対称であるべきです。 もし設定した値を tclean があまりよくないと判断したならば、ログウィンドウでもっと良い値を提示するでしょう。
tclean タスクは対話的モードで開始され、閾値、メジャーサイクル、マスクを手動でコントロールできます。 以前の clean に慣れ親しんだユーザーにとっては、clean タスクとは対照的に、niter (最大反復回数) を設定しないままにしておくと、 niter = 0 がデフォルトとなり、tclean はクリーニングをおこなわないことに注意してください。 モデル内にクリーンなコンポーネントがない場合、例えば、与えられた MS に対して、このフィールドとスペクトルウィンドウに対して tclean を最初に起動した場合、作成されるのはダーティキューブまたはイメージです。ここでは、まず niter を設定しない状態で tclean を実行してみます。

# In CASA
tclean(vis='sis14_twhya_calibrated_flagged.ms',
       imagename='phase_cal',
       field='3',
       spw='',
       specmode='mfs',
       deconvolver='hogbom',
       gridder='standard',
       imsize=[128,128],
       cell=['0.1arcsec'],
       weighting='briggs',
       threshold='0.0mJy',
       interactive=True)

CASA は次のような、クリーニングをしていないことを示す警告を出すはずです。


WARN task_tclean	Restoring with an empty model image. Only residuals will be processed to form the output restored image.

もう一度コマンドを実行し、今度は niter を大きな数字に設定します。

# In CASA
os.system('rm -rf phase_cal.*')

tclean(vis='sis14_twhya_calibrated_flagged.ms',
       imagename='phase_cal',
       field='3',
       spw='',
       specmode='mfs',
       deconvolver='hogbom',
       gridder='standard',
       imsize=[128,128],
       cell=['0.1arcsec'],
       weighting='briggs',
       threshold='0mJy',
       niter=5000,
       interactive=True)

図 3: tclean の GUI では、クリーンマスク (白い楕円) が明るい領域を囲んでいます。

tclean タスクを実行すると、図 2 のような GUI が表示されます。 tclean のビューワーで、ボタンが新しい楕円形のマスク領域を追加するように設定されていることを確認します。 楕円の中に "R" が表示されているアイコンをクリックする必要があるかもしれません。 点光源である二次較正天体が選択されているので、この時点でフィールドの中央に点光源が見えるはずです。 シグナルがある部分 (中央の点だけ) の周りに楕円形のマスクを描きます。楕円の内側をダブルクリックして、白くなるのを見ます (図 3 参照)。 マスクを設定するときは、シグナルの部分のみを選択し、それ以外はあまり含まないようにします。GUI の "Next Action" セクションに、幾つかのコントロールボタンがあります。 緑色の丸いボタンを押すと、クリーンが開始されます。これでクリーンのサイクルが実行され、その後戻ります。最初の clean の後、残差が表示されます。 残差を光源と比較し、十分だと判断したら(または閾値 (デフォルトでは 0 [mJy]) を満たしたら、つまり残差が負になったら)、赤い "X" をクリックすると tclean は終了します。 この例では、2 回のクリーンでうまくいくはずです。 より複雑なターゲットの場合何回もクリーンを行う必要があり、各サイクルの後に、残差の見た目に基づいて新しいマスク領域を追加することが可能です。 tclean でできたファイルを確認します。

# In CASA
ls

"phase_cal" で始まる、拡張子だけが違ったファイルが幾つかできています。
  1. .image : tclean した最終的なイメージ
  2. .mask : (最後のクリーンのサイクルで) 使用したマスクのイメージ
  3. .model : tclean で用いた (抽出した) clean components (単位は [Jy/pixel])
  4. .pb : 主ビームのレスポンス
  5. .residual : 最終的な残差マップ (最終イメージの "dirty" 部分に相当)
  6. .psf : 合成ビームのイメージ (dirty beam)
  7. .sumwt : 重みの和を含む単一ピクセルのイメージ
CASA ビューワーを使うとこれらを見ることができます。ビューワーの起動は viewer() 又は、imview() でおこないます。 UNIX / Linux のコマンドプロンプト (CASA 外) からは casaviewer でスタンドアローンのビューワーを起動することもできます。

# In CASA
imview("phase_cal.image")

ビューワーを使って対話的に他の画像も読み込んでみましょう。 ビューワーにある "Images" 用のテープデッキ・コントローラ (進む、戻る、止まるなどのボタン) を使って変更できます。

tclean のパラメーターを変更する

tclean にはたくさんのオプションがあります。inp tcleanと入力すれば、タスクの入力リストを見ることができます。 よく変更されるオプションの 1 つは、UV データをフーリエ平面画像に展開するのに使われる重み付けのパラメーターです。 この重み付けは、最初の例では robust=0.5 でした (デフォルト)。robust パラメータの値を -2 (uniform) と 2 (natural) の間でいくつか変えてみてください。 ビームサイズがどのように変化するかと、また最終画像のノイズに注目して下さい。 ノイズレベルは、ビューワーを使ってチェックできます。 画像上にボックスを描き、"Regions" ボックスの "Statistics" タブでその領域の画像の統計値を見てください ('Regions' ボックスが表示されない場合は、ビューワー GUI 上部の "View" ドロップダウンメニューからオンにしてください)。 以前に実行したことがある場合は、古いバージョンの画像を削除してください。

# In CASA
os.system('rm -rf phase_cal_robust.*')

briggs weighting と robust = -1 を指定して tclean を呼び出します。 クリーンマスクを作成し、残差のレベルに満足するまで、クリーンを数サイクル実行します。

# In CASA
tclean(vis='sis14_twhya_calibrated_flagged.ms',
       imagename='phase_cal_robust',
       field='3',
       spw='',
       specmode='mfs',
       gridder='standard',
       deconvolver='hogbom',
       imsize=[128,128],
       cell=['0.1arcsec'],
       weighting='briggs',
       robust=-1.0,
       threshold='0mJy',
       niter=5000,
       interactive=True)

図 4: Briggs による重み付けパラメーターを設定しての tclean をおこなった位相較正天体のビューワー画像。

結果を確認します。

# In CASA
imview("phase_cal_robust.image")

他の較正天体 (フィールド 0 と 2 。listobs 出力で "ScanIntent" をもう一度確認してください) をイメージングして、画像サイズとセルサイズを大きくしたり小さくしたりしてみてください。 例えば、わずかに分解された一次 (フラックス) 較正天体である Ceres (フィールド 2; "ScanIntent" = "CALIBRATE_AMPLI#ON_SOURCE") を見てみましょう:

# In CASA
os.system('rm -rf amp_cal_robust.*')

tclean(vis='sis14_twhya_calibrated_flagged.ms',
       imagename='amp_cal_robust',
       field='2',
       spw='',
       specmode='mfs',
       gridder='standard',
       deconvolver='hogbom',
       imsize=[128,128],
       cell=['0.1arcsec'],
       weighting='briggs',
       threshold='0mJy',
       niter=5000,
       interactive=True)

imview("amp_cal_robust.image")

上の plotms で見たように、ケレスはある程度分解されています。 ピクセルサイズを非常に大きくしてみると、画像が壊れるのがわかるでしょう。 tclean (tclean はデコンボリューションをピクセル単位で量子化する) の目的から、ピクセルサイズは合成ビームに比べて小さくすることをお勧めします。 ピクセルサイズが合成ビームに比べて大きい場合、tclean の方法とは関係なく、一般的に画像は劣化します。

# In CASA
os.system('rm -rf amp_cal_bigpix.*')

tclean(vis='sis14_twhya_calibrated_flagged.ms',
       imagename='amp_cal_bigpix',
       field='2',
       spw='',
       specmode='mfs', 
       gridder='standard',
       deconvolver='hogbom',
       imsize=[32,32],
       cell=['0.5arcsec'],
       weighting='briggs',
       threshold='0mJy',
       niter=5000,
       interactive=True)

imview("amp_cal_bigpix.image")

図 5: 非常に大きなピクセルサイズを設定して tclean をおこなった場合の振幅較正天体。

設定したピクセルサイズと最終的なビームサイズを比較して、どのような画像になるかを確認してみましょう。

較正とフラグの効果

ここでは、各種較正とフラグがどのような効果を与えるかを説明します。 イメージングの工程だけを見たい場合にはスキップできます。

ここで扱っているデータは、様々な較正過程を経ています。これらの較正はどのような効果があったのでしょうか? 較正をおこなった場合とおこなわなかった場合、フラグを立てた場合と立てなかった場合で、二次較正天体の最終画像がどのように変わってくるかを試してみます。 まず、未較正のデータを作業ディレクトリーからコピーする必要があります。 これはデータを保存している場所によって、どこからコピーしてくるかが変わりますが、最終的に未較正のデータが現在の作業ディレクトリーにあれば問題ありません。 完全なデータパッケージをダウンロードしている場合は、以下のコマンドによって未較正データを作業ディレクトリーにコピーします。

# In CASA
os.system("tar -xvf ../working/sis14_twhya_uncalibrated.ms.tar")
os.system("rm -rf sis14_twhya_uncalibrated.ms")
os.system("cp -r ../working/sis14_twhya_uncalibrated.ms .")

較正済みデータセットのみをダウンロードした場合は、以下のように wget を使用して未較正データをダウンロードすることができる。

# In bash
% wget -r -np -nH --cut-dirs=4 --reject "index.html*" https://bulk.cv.nrao.edu/almadata/public/working/sis14_twhya_uncalibrated.ms/

又は、このリンクからブラウザ経由でファイルをダウンロードすることもできます。

図 6: 較正を何もしていない位相較正天体のビューワー画像

未較正データを tclean してみます。 二次較正天体 (フィールド 3) に対して、前と同じコマンドを実行してみます。

# In CASA
os.system('rm -rf phase_cal_uncalibrated.*')
tclean(vis='sis14_twhya_uncalibrated.ms',
       imagename='phase_cal_uncalibrated',
       field='3',
       spw='',
       specmode='mfs',
       gridder='standard',
       deconvolver='hogbom',
       imsize=[128,128],
       cell=['0.1arcsec'],
       weighting='natural',
       threshold='0mJy',
       niter=5000,
       interactive=True)

クリーンマスクをかけようとする天体が見つかればよいのですが…。

# In CASA
imview("phase_cal_uncalibrated.image")

生データ (それでも Tsys と WVR は補正されています)では、フィールド全体に較正天体のエコーが見えます。

次はフラグありなしによる違いです。フラグのついていないデータを作業ディレクトリーからコピーしてきます。

# In CASA
os.system("tar -xvf ../working/sis14_twhya_calibrated.ms.tar")
os.system("rm -rf sis14_twhya_calibrated.ms")
os.system("cp -r ../working/sis14_twhya_calibrated.ms .")

wget でダウンロードすることもできます。

# In bash
% wget -r -np -nH --cut-dirs=4 --reject "index.html*" https://bulk.cv.nrao.edu/almadata/public/working/sis14_twhya_calibrated.ms/

ブラウザ経由ならば、このリンクからダウンロードできます。 前と同じパラメーターで二次較正天体の、較正済み、かつ、フラグ処理をしていないデータを tclean してみます。

# In CASA
os.system('rm -rf phase_cal_unflagged.*')

tclean(vis='sis14_twhya_calibrated.ms',
       imagename='phase_cal_unflagged',
       field='3',
       spw='',
       specmode='mfs',
       gridder='standard',
       deconvolver='hogbom',
       imsize=[128,128],
       cell=['0.1arcsec'],
       weighting='natural',
       threshold='0mJy',
       niter=5000,
       interactive=True)

imview("phase_cal_unflagged.image")

図 7: 悪いデータのフラグ処理をしていない位相較正天体のビューワー画像

較正されていないデータとは違って、フラグ処理をしていないデータは幾分まともですが、残差に明らかに人工的な模様が見えます。 フラグ処理をすることによって最終的なイメージの質は向上しますが、今回のような全体に質のよいデータでは、フラグ処理をしていなくてもターゲットの目的天体も一応見えています。

目的天体のイメージング

様々な較正の過程は、目的天体のデータ (サイエンスデータ) を較正することにあります。ここでは、TW Hydra の観測はフィールド 5 がサイエンスデータになります。 基本的なイメージングのチュートリアルの最後のステップは、TW Hydra の連続波イメージングです。 まず、CASA のタスク split を使って、サイエンスデータをデータセットに "分割" します。 これは必ずしも必要ではないですが、これはデータの管理を容易にする一般的なステップです。 同時に、width=8 を使用してデータを平滑化 (平均化) します。この時点では連続波イメージングが対象のため、多くの情報を失うことなくデータ量を減らすことができます。 このデータセットはすでに管理しやすいように設計されていますが、ALMA は非常に大きなデータセットを生成することがあるので、このスムージングのトリックは役に立つでしょう。

# In CASA
os.system('rm -rf twhya_smoothed.ms')

split(vis='sis14_twhya_calibrated_flagged.ms', field='5', width='8', outputvis='twhya_smoothed.ms', datacolumn='data')

listobs('twhya_smoothed.ms')

次に、分割したデータの連続波イメージを作成します。 TW Hydra は、新しいデータセットではフィールド 0 として再ラベル化されています。 ここでも多周波合成モード ("mfs") を使い、ピクセルサイズを上よりもやや小さく (=イメージサイズをやや大きく) します。 これは TW Hydra は広がっている天体であるため、そして、"briggs" 重み付けを使うためビームがやや小さくなるためです。 ここでは対話的モードを指定し、閾値はとりあえず設定しないことにします。

# In CASA
os.system('rm -rf twhya_cont.*')

tclean(vis='twhya_smoothed.ms',
       imagename='twhya_cont',
       field='0',
       spw='',
       specmode='mfs',
       gridder='standard',
       deconvolver='hogbom',
       imsize=[250,250],
       cell=['0.08arcsec'],
       weighting='briggs',
       robust=0.5,
       threshold='0mJy',
       niter=5000,
       interactive=True)

図 8: TW Hydra の連続光放射

マスクツールを使って天体の周りにクリーンマスクを描き、TW Hydra ディスクからの放射がその周りの残差以下か同等になるまでクリーンを実行します。 この対話的なクリーニングモードでは、いつクリーンを停止させるかを決め、その時点で赤い X のボタンを押します。 代わりに、閾値まで自動的にクリーニングしたい場合は、上記の tclean の呼び出しの際に、或いはビューワーウィンドウに入力して、RMS ノイズの数倍の閾値を指定できます。 画像を見てください。TW Hydra は非常に明るく、ビームの大きさに対して大きいのがわかります。残差は完璧ではありませんが、この後の過程で改善していきます。

# In CASA
imview("twhya_cont.image")

しかし、待ってください。 このスペクトルウィンドウにはスペクトル線 (輝線) があると言いませんでしたか?

はい、今回連続光として撮像したスペクトルウィンドウには、N2H+ の輝線も含まれています。 この場合、N2H+ の輝線は暗いので、イメージング前に輝線のチャンネルにフラグを立てなくても、最終的な連続体画像にほとんど違いはありません。 このチュートリアルと他の連続光ファーストルック・チュートリアルでは、イメージングとセルフキャリブレーションの基本的なステップに重点をおいているため、輝線を無視しています。 しかし、一般的には、連続体をイメージングする前に、自分のデータで輝線を含むチャンネルにフラグを立てるべきです。 スペクトルウィンドウで連続光をイメージングする前に、輝線にフラグを立てる方法については、 IRAS16293_Band9_-_Imaging (英語) の上級チュートリアルを参照してください。

非対話的モードの tclean

tclean は、対話的なガイダンスなしで実行するように実行することもできます。 指定する 3 つの主なパラメーターは、停止する閾値 (クリーンする領域の最大残差がこの閾値より低いとき、クリーンのサイクルがは停止します)、 マスク (クリーンの工程がシグナルを識別する領域)、および最大反復回数です。 最大反復回数は厳密には必要ではなく、一般的にはフェイルセーフとして使用することが推奨されます。 つまり、もしクリーンのサイクルがその回数に到達したら何かが間違っている、というくらいの数値に設定するのがよいでしょう。

まず、クリーンマスクを決定する方法を試してみましょう。 先ほど作った画像を見てください。天体からのものと思われる明らかな放射は全て、(100,100) から (150,150) のような範囲で囲まれたボックスに含まれています。 tclean の呼び出しで "mask "パラメータを使って、このボックスをマスクに設定します。 また、ファイル (たとえば、以前の対話的バージョンの tclean で作成されたもの) を指定して設定することもできます。

次に、tclean の停止閾値を決めましょう。 もう一度、imview を使って前の画像を開き、ノイズを推定するために、天体からかなり離れたところにボックスを描きます。 先ほどと同じように、ビューワーの "Regions" ボックスの "Statistics" タブを見て、RMS ノイズを求めます。 7 [mJy/beam] のようなものが見えます。閾値をこの約 2 倍、~15 [mJy/beam] に設定します。 通常、デコンボリューションされた画像に偽の天体が追加されるのを避けるために、rms ノイズの数倍の閾値が推奨されます。 つまり、tclean がランダムなノイズスパイクを天体として扱い、画像からそれをデコンボリューションすることを避けたいからです。 これは、セルフキャリブレーションを行っている場合に特に問題となります。 最後に niter=10000 を設定します。この回数に達する前にクリーンが終了することを期待します。 これはクリーンが永遠に実行されないようにするための単なる大きな数字です。

# In CASA
os.system('rm -rf twhya_cont_auto.*')

tclean(vis='twhya_smoothed.ms',
       imagename='twhya_cont_auto',
       field='0',
       spw='',
       specmode='mfs',
       gridder='standard',
       deconvolver='hogbom',
       imsize=[250,250],
       cell=['0.08arcsec'],
       mask='box [ [ 100pix , 100pix] , [150pix, 150pix ] ]',
       weighting='briggs',
       robust=0.5,
       threshold='15mJy',
       niter=10000,
       interactive=False)

imview('twhya_cont_auto.image')

この非対話的モードは、多くの時間を節約でき、非常に再現性が高いという利点があります。 interactive=True でクリーンプロセスを開始し、右上の青い矢印ボタンをクリックすると、"ハイブリッド" モードも利用できます。 これは tclean に最大反復回数または閾値に達するまで処理を進めるように指示する。 この組み合わせモードは、tclean に使用するマスクを手動で描画できるところが利点です。 ビューワー GUI で、閾値と最大反復回数 (メジャーサイクル数とサイクルあたりの反復回数の積) の両方を手動で設定することもできます。 ただし、較正が不確かな画像、特に光源が明るい画像については、少なくとも初回は対話的にクリーニングするのがよいでしょう。 画像が "ダイナミックレンジに制限がある" (すなわち、品質が較正とデコンボリューションの精度によって設定される) 場合、正しいしきい値を予測することは困難です。 クリーニングのサイクルが進むにつれて、画像の残差は減少するはずです。 ロガーの tclean に関する部分でこれを確認できます。残差が増加し始めたら、停止しきい値を変更する必要があるかもしれません。

プライマリービーム補正

tclean の重要な点は、デフォルトでは tclean によって生成されたイメージが、アレイ内の個々のディッシュのプライマリービーム (視野) に対して補正されないことです。 プライマリービームの応答は、通常、視野の中心で値 1 のガウスです。 天文学的に正しい空の画像を形成するには、tclean の出力をこのプライマリービーム (モザイクの場合は、モザイクを作成するために使用されるプライマリービームパターンの組み合わせ) で分割する必要があります。tclean の場合、これを行うには 2 つの方法があります。
まず、tclean の実行時にパラメータ pbcor = True を設定する方法です。 これにより、.image.pbcor という拡張子を持つ追加画像が生成され、これがプライマリービーム用に補正されたクリーニング画像となります。
もう一つの方法は、CASA はこの補正に必要なプライマリービーム情報が保存された拡張子 .pb のファイルを用いる方法です。 CASA のタスク impbcor を使用して、.pb 画像と tclean からの出力画像を結合し、プライマリービーム補正された画像を生成することができます。
まず、古いプライマリービーム補正画像があれば、それを削除します。

# In CASA
os.system('rm -rf twhya_cont.pbcor.image')

次にイメージを補正します。

# In CASA
impbcor(imagename='twhya_cont.image',
        pbimage='twhya_cont.pb',
        outfile='twhya_cont.pbcor.image')

出力されたイメージをチェックします。

# In CASA
imview('twhya_cont.pbcor.image')

プライマリービーム補正前のイメージで作業することは、ノイズがフィールド全体で同じであるため、非常に便利です (例えば、これはシグナルを探索するためのクリーンなデータセットです)。 しかし、サイエンスのためにフラックスや強度を計算する前に、この補正を適用することを忘れないことが非常に重要です。

図 9: プライマリービーム補正された、TW Hydra の連続光イメージ



ページ先頭に戻る
CASA guides "First Look at Imaging CASA 6.4" へ (本家のサイト:英語)

Last Update: 2023.09.30