A-4.ビジビリティー振幅の較正

AIPSによる処理(4) 振幅の較正

いよいよデータの較正にとりかかります。較正とは、観測された生データから真の入力値を推定するために、伝達関数のパラメーターを求めることです。一般に、私たちが手にする観測量は、観測装置に入力された値が伝達関数の影響を受けて出力されたものです。入力と出力との関係のことを伝達関数と呼びます。詳しくは干渉計サマースクール2005のこのページをご覧下さい。
たとえば干渉計で測定されるビジビリティは、真のビジビリティにアンテナ複素ゲインを掛け合わしたものです。式で書くと、

となります。ここで、tは時刻,νは周波数, i, jはアンテナ番号, gはアンテナの複素ゲイン, Vは真のビジビリティ,

は観測されたビジビリティです。
正しいビジビリティを推定するには、アンテナ複素ゲインの補正が必要です。アンテナ複素ゲインはその名の通り複素数で、振幅と位相の項があります。AIPSでは振幅の較正と位相の較正を別々に行います。どちらを先にやってもよいのですが、ここではまず振幅の較正を行います。つまり、アンテナ複素ゲインの振幅項を求めよう、というわけです。

1. アンテナゲインの振幅項

相関器から出力されるビジビリティは正規化された相関関数をフーリエ変換したものですので、SEFD (System Equivalent Flux Density : システム雑音等価フラックス密度単位=Jy) を掛けるとFlux Density (フラックス密度: 単位 = Jy)になります。つまり、

です（平方根をとるのは、二つのアンテナゲインの幾何平均が基線の合成SEFDになるから）。
SEFDは、システム雑音温度Tsysとアンテナ有効開口面積Aeを使って

と表されます。ここで、 $k_{B}$ = 1.38× $10^{3}$ はボルツマン定数（Jy= $10^{-26}$ WHz $^{-1}$ $m^{-2}$ なので、それに合わせてボルツマン定数も $10^{26}$ 倍している）です。定数2が付くのは、偏波の1成分だけを受信している状況を考えているからです。結局、ゲインの振幅項を得るにはアンテナ毎に有効開口面積Aeとシステム雑音温度Tsysの情報があればよい、というわけです。
Tsysは大気込みの雑音温度で、天候に左右されますから、時間の関数です。周波数によってもTsysは変化しますけど、IF内の狭い帯域の中ではほぼ一定と思っても構わないので、AIPSではTsysを時間だけの関数として扱います。観測時に測定されたTsysデータは、TY extension tableに記録してあります。
有効開口面積は時間に依存するというよりは、アンテナの仰角に依存します。重力変形が有効開口面積減少をもたらすからです。そこで、有効開口面積をAIPSの中では仰角 (EL) の多項式関数として扱い、多項式の係数を GC extension tableに記録してあります。なお、AIPSではAeの値ではなく、Ae/2kBの値をGC extension tableで扱います。

2. AIPSにおける較正のポリシー

2.1 データは神聖なり

AIPSにおける較正のプロセスは、「データは神聖なり」というポリシーを貫いているという特徴があります。これは、観測されたオリジナルのビジビリティセットには一切手を加えずに、そのまま保存し続ける、という方針のことです。それでは較正されたビジビリティはどのように作成するのでしょうか。AIPSにおける較正とは、較正テーブルを作成すること、にほかなりません。
較正テーブルにはアンテナ複素ゲイン $g^{\frac{-1}{2}}$ が記録されています。この $g_{i}^{\frac{-1}{2}}$ , $g_{j}^{\frac{-1}{2}}$ を観測されたビジビリティ

i, jに掛けることで、較正されたビジビリティ $V_{i,j}$ が表示される（あるいは出力される）というわけです。観測されたビジビリティデータは一切変更されることなく、保存されます。

AIPSにおけるビジビリティデータとCLテーブルの扱い→

この方式の利点は、

大容量になりがちなビジビリティデータを書き換えたり複製する必要がなくなるので、計算負荷や記憶容量を節約できます。
較正テーブルを気軽に試行錯誤できます。較正テーブルの作成に失敗することはよくあることです。チェックのために表示して、「あ、まちがった」と思ったら、較正テーブルを新たに作り直せばよいわけです。複数の較正テーブルを保存しておいて、どのバージョンの較正テーブルを最終的に使用するかを選択できます。
較正テーブルの作成をstep-by-stepにできます。第1段階では正規化, 第2段階で振幅の時間変化較正, 第3段階で位相の時間変化の較正, 第4段階で帯域通過特性の較正…と、較正テーブルを順次アップデートしていくのが容易にできます。

2.2 CLテーブルとBPテーブル

「較正テーブル」と呼んだ複素ゲインg(ν, t)を、AIPSでは周波数に依存する項と時間変化する項に分離して扱います。

というように、B(ν)を帯域通過特性としてBP extension tableに、G(t)をCL extension tableに、それぞれ格納します。AIPSにおいてBP tableは時間変化してもよいのですが、観測中に帯域通過特性は一定として扱ってもたいていは大丈夫なので、このコースでは時間不変なものとして扱います。下図に、BP table と CL table の例を示します。

BP tableの例。横軸が周波数分光チャネル番号で、紫は位相、緑は振幅を表している。各アンテナ毎に4つのIFがある観測なので、その4IF分を表示している。IF帯域の高周波端で振幅が低下しているのはフィルターの特性を表している。

CL tableの例。横軸が時間で、縦軸に振幅項をプロットしている。10個のアンテナについて、IF=1におけるゲインを表示している。縦軸が約940 milligain になっているのは、値が0.94程度ということ。この段階ではSEFDはまだ反映させておらず、自己相関関数を用いた正規化だけを行っているので、量子化損失の値だけが反映されて1に近い値になっている。

2.3 較正のプロセス

それでは、このコースにおける具体的な較正の過程を俯瞰してみましょう。下図をご覧下さい。

較正のプロセス。CL.*とあるファイルはCL extension tableで、較正の段階を経るにつれてバージョンが1から4まで進む。GCはgain curve extension tableで、Ae/2kBの値をELの関数として格納してある。TYはTsys extension tableでシステム雑音温度の値が時間の関数として記録されている。SN.* はSolution extension tableで、観測されたビジビリティを元に解を得て作成したり、GCやTYのデータを元に作成する。このSNをCLに適用することで、CL tableを更新することができる。BPはBandpass extension tableで、帯域通過特性を格納する。ACCOR, CLCAL, APCAL, FRING, BPASS, SPLITはそれぞれ使用するtaskの名前である。

上図のように、較正のプロセスはCL tableを更新してバージョンアップしていくことです。CLのバージョン1は、「何もしない」すなわち1という値がひたすら記録されています。1をビジビリティに掛け合わせても値は変化しませんからね。バージョン1のCL talbleからスタートして、正規化係数を適用したり、SEFDを適用したり…と進み、最終的なバージョンのCL tableを適用すると、較正されたビジビリティが得られる、というわけです。
CLの更新に用いられるのがSN tableで、これはビジビリティに対してある解を求めるような方程式を解いたり、GCやTYを元に作成したりします。
CLは観測時間内で一定の時間刻みでまんべんなく較正データが記録されているのに対して、SNには不連続的に解が記録されます。そのため、SNをCLに適用するCLCALというtaskにおいては、補間や平滑化という操作が必要になります。

3. 自己相関データを用いて正規化する (ACCOR)

観測されたビジビリティは、相関器の中であるていどの正規化はされていますが、自己相関データを用いることで厳密な正規化がなされます。この操作をaccor という task で行います。
観測されたビジビリティ

#ref error ：ご指定のファイルが見つかりません。ファイル名を確認して、再度指定してください。 (v.png)

i, jは、相互相関関数 $C_{i,j}(\tau)$ をフーリエ変換して得られます。

相互相関関数の正規化は、自己相関関数のτ=0における値の幾何平均で割り算することです。

従って、 $C_{ii}^{\frac{-1}{2}}(0)$ の値を較正テーブルに記録することで、ビジビリティを正規化できるわけです。accorは $C_{ii}^{\frac{-1}{2}}(0)$ の平均値を与えられた時間きざみ毎に算出して、SN extension tableに出力します。
それでは、実際にaccorを実行してみましょう。

task 'accor'ACCORというtaskの使用宣言
getn 2カタログ番号2番のファイルを選択
AIPS 1: Got(1) disk= 1 user=3018 type=UV >BK084.MSORT.1
timer 0全ての時間範囲で正規化を行う
solin 1　　時間きざみ幅を1分に設定
inp　　パラメーターの一覧を表示して確認します。

（accorのパラメーターの一覧はこちら）
goと打って実行します。終了したら、imhで SN extension tableのバージョン1が追加されたことを確認しましょう。

3.1 できたSN table を確認する (SNPLT)

SN extension table をTV画面に表示して確認するには、SNPLT という task を使用します。

task 'snplt'　SNPLTというtaskの使用宣言
getn 2　　カタログ番号2番のファイルを選択
AIPS 1: Got(1) disk= 1 user=3018 type=UV >BK084.MSORT.1
inext 'sn'　　SN extension tableを表示する、と指定 (SNPLTはCLやTYも表示できる)
inv 1SN extension tableのバージョン1を表示するように指定
optyp 'amp'振幅を表示するように指定
nplot 101ページの画面に10個（アンテナの数）のデータを表示す>るよう指定。
inp　　パラメーターの一覧を表示して確認します
（snpltのパラメーターの一覧はこちら）
tvinit　　TV画面を初期化（クリア）

goと打って実行すると、TV画面上に下図のようなプロットが表示されます。横軸が時刻で、縦軸が振幅です。milligain 単位で表示されていることに注意… 例えば960 milligain = 0.96で、ほぼ1に近い値となります。

3.2 SN table をCL tableに適用する (CLCAL)

SN extension table に格納されたゲインをCL extension tableに適用します。それにはCLCAL という task を使用します。
CL tableにはゲインが0.1分の時間きざみ幅 (この時間きざみ幅はINDXRを実行したときに指定しました)ですき間なく記録されているのに対して、SN tableに格納された解は1分間隔です。従って、CLCALを実行するに当たっては補間の処理が必要になります。補間の方法について、CLCALの入力パラメーターで指定してやる必要があります。

task 'clcal'CLCALというtaskの使用宣言
getn 2　　カタログ番号2番のファイルを選択
AIPS 1: Got(1) disk= 1 user=3018 type=UV BK084.MSORT.1
source 'da193'''　　 DA 193を観測していた時間帯のCLテーブルデータだけに適用
calsou 'da193'''　　DA 193を観測していた時間帯のSNテーブルデータだけを使用
freqid 1　　周波数番号を1番 (15.4 GHz) に指定
INTERPOL 'self'　　天体番号ごとそれぞれに補間
SAMPTYPE 'mwf'　　補間方法として median window filterを使用
BPARM 1　　median window filterの時間窓を1 hourに
SMOTYPE 'ampl'振幅の平滑化
SNVER 1　　SN extension tableのversion 1を使用
GAINVER 1　　CL extension tableの version 1を補正元に指定
gainu 2　　補正後の結果をCL extension table version 2に格納
inp　　パラメーターの一覧を表示して確認します。

（clcalのパラメーターの一覧はこちら）
goと打って実行します。終了したら、imhで CL extension tableのバージョン2が追加されたことを確認しましょう。

3.3 できたCL table を確認する (SNPLT)

CL tableもSN tableと同様にSNPLTというtaskで表示・確認できます。

tget snpl　　SNPLTの使用宣言。tget という verbを用いると、前回実行時のパラメーターがセットされる。
inext 'cl'　　CL extension tableを表示するよう指定
inv 2　　CL extension tableのバージョン2を表示するように指定
optyp 'amp'　　振幅を表示するように指定
inp　　パラメーターの一覧を表示して確認します。
（CL確認のためのsnpltのパラメーター一覧はこちら）
tvinit　　TV画面を初期化（クリア）

goと打って実行すると、下図のようにTV画面に表示されます。

4. SEFDの補正をする (VLBAMCAL, SETJY, APCAL)

正規化されたビジビリティにSEFDの幾何平均を掛けると、Jy単位のビジビリティに較正されるのでした。SEFDはアンテナゲインとシステム雑音温度 (Tsys) から求めるのでしたよね。アンテナゲインはGC extension tableに、TsysはTY extension tableにそれぞれ記録されています。これらの情報を元にSEFDを計算してSN tableに格納するには、3つの手順を踏みます。VLBAMCAL というマクロ、それにSETJYおよびAPCALという tasks です。

4.1 VLBAMCALでTYとGCを準備

VLBAMCALまマクロです。マクロとは複数のtaskやverbをひとまとめに記録したスクリプトで、実行するとスクリプトに記述された順番にtaskやverbを次々と実行してくれます。VLBAMCALは、GC extension tableとTY extension tableの調整をしてくれるマクロです。

run vlbautil　　VLBAMCALというマクロを実行するために、VLBAUTILというライブラリーを起動
getn 2　　カタログ番号2番のファイルを選択
AIPS 1: Got(1) disk= 1 user=3018 type=UV BK084.MSORT.1
inp vlbamcal　　VLBAMCALの入力パラメーターを確認
AIPS 1: VLBAMCAL: Procedure to merge redundant calibration data
AIPS 1: Adverbs Values Comments
AIPS 1: ----------------------------------------------------------------
AIPS 1: INNAME 'BK084' Input file name
AIPS 1: INCLASS 'MSORT' Input file class
AIPS 1: INSEQ 1 Input file sequence number
AIPS 1: INDISK 1 Disk number for input file
AIPS 1: BADDISK *all 0 List of disks not to be used
AIPS 1: for scratch files.
AIPS 1:
AIPS 1: VLBAMCAL is defined in the VLBAUTIL run file.
vlbamcal　　VLBAMCALを実行

たくさん現れるメッセージは「何をやっているか」を表すものですが、特に気にしなくてもOKです。終了した時に、GCとTYそれぞれのバージョン1があることを、imhで確認してください。

4.2 SETJYでシステム雑音温度測定時の天体のフラックス密度を指定

TY extension tableにはシステム雑音温度 (Tsys) が記録されています。VLBAでは、天体の信号を受信している観測中に並行してTsysを測定するので、実は測定されているのはTsysでなく Tsys + Ta です。ここで、Taは天体のアンテナ温度で、フラックス密度Sの天体を観測した時のアンテナ温度は Ta = Ae S / $2k_B$ という関係です。あまり強くない天体では Tsys ≫ Ta となるのでTaは無視してもいいのですが、厳密に較正を行うためには天体のフラックス密度をAIPSに教えてあげて、Taの補正をする必要があります。AIPSでは天体のフラックス密度はSU extension tableに記載されます。記載するにはSETJYというtaskを使用します。
このコースではDA 193という天体を観測していますので、DA 193のフラックス密度をSETJYで記録しましょう。さて、DA 193のフラックス密度はいくらでしょう？NRAO (米国国立電波天文台）の較正天体モニター観測や、UMRAO（ミシガン州立大)のモニター観測データが便利です。ここではUMRAOのデータを用いましょう。それによると、15 GHzにおいてDA 193のフラックス密度は4.8 Jyです。

task 'setjy'　　SETJYというtaskの使用宣言
getn 2　　カタログ番号2番のファイルを選択
AIPS 1: Got(1) disk= 1 user=3018 type=UV BK084.MSORT.1
source 'da193'''　　天体名をDA193と設定
zerosp 4.8, 0　　フラックス密度を4.8 Jyと設定
freqid 1　　周波数番号を1番 (15.4 GHz) に指定
optyp ''　　OPTYPEのパラメーターをクリア
inp　　パラメーターの一覧を表示して確認します。

（setjyのパラメーターの一覧はこちら）
goと打って実行すると、メッセージ画面にこのようなメッセージが現れて終了します。

4.3 APCALでSEFDをSN extension tableに書き込み

やっとTY, GC, SUが揃ったのでSEFDを計算できるようになりました。APCALというtaskでSEFDを計算してSN extension tableに書き込みます。

task 'apcal'　　APCALというtaskの使用宣言
getn 2　　カタログ番号2番のファイルを選択
AIPS 1: Got(1) disk= 1 user=3018 type=UV BK084.MSORT.1
source 'da193'''　　天体名をDA193と設定
freqid 1　　周波数番号を1番 (15.4 GHz) に指定
opcode 'LESQ'　　最小二乗法でフィット
dotv 1　　結果をTV画面に表示
inv 1　　WX extension table (気象データ) のversion 1を使用
tyv 1　　TY extension table (Tsysデータ) のversion 1を使用
gcv 1　　GC extension table (アンテナゲインデータ) のversion 1を使用
snv 2　　計算結果をSN extension table のversion 2に書き込み
dofit 1, 0　　大気の光学的厚みの補正を行う
inp　　パラメーターの一覧を表示して確認します。

（apcalのパラメーターの一覧はこちら）
goと打って実行します。終了したらimhで SN extension tableのバージョン2が追加されたことを確認しましょう。

APCALがうまくいかなかった場合の対処: ANTABを使用

AIPSのバージョンが31DEC06より古い場合、あるいは使用するデータの相関処理が古い場合、VLBAMCALで作成したGCテーブルをAPCALでうまく使えない場合があるようです。その場合には、ANTABというtaskを利用して、テキストファイルからGCテーブルを作成し、APCALを実行します。
まず、下記URLからアンテナゲイン情報が書き込まれたテキストファイル BK084.GAIN をダウンロードして、AIPS_ROOT/FITS に保存してください。AIPS_ROOTが /usr/aips であれば、/usr/aips/FITS/BK084.GAIN というファイルに保存するわけです。

http://astro.sci.kagoshima-u.ac.jp/omodaka-nishio/member/kameno/AIPS/BK084.GAIN : 739-byte ASCII file

次にANTABというtaskで、テキストファイルからGCテーブルを作成します。

task 'antab'　　APCALというtaskの使用宣言
getn 2　　カタログ番号2番のファイルを選択
AIPS 1: Got(1) disk= 1 user=3018 type=UV BK084.MSORT.1
infile 'FITS:BK084.GAIN'　　読み込むテキストファイル名を指定
gcv 2　　AIPSで作成されるGC extension tableのバージョン
tyv 2　　Tsys情報を含まないテキストファイルを読み込むので、このTYはダミー
inp　　パラメーターの一覧を表示して確認します。

goと打って実行します。終了したらimhで GC extension tableのバージョン2が追加されたことを確認しましょう。
これでAPCALが実行できるようになりました。

tget apcal　　さっきAPCALを実行したときのパラメーターを再利用するため、tgetします
gcv 2　　新たに作成した、version 2のGC extension tableを使います

goと打って実行します。終了したらimhで SN extension tableのバージョン2が追加されたことを確認しましょう。
以上で ANTAB による回避は終わりです。では、コースに復帰しましょう。

4.4 APCALでできたSN table を確認する (SNPLT)

作成されたSN version 2をSNPLT で表示して確認しましょう。

task 'snplt'　　SNPLTというtaskの使用宣言
getn 2　　カタログ番号2番のファイルを選択
AIPS 1: Got(1) disk= 1 user=3018 type=UV BK084.MSORT.1
inext 'sn'　　SN extension tableを表示する、と指定 (SNPLTはCLやTYも表示できる)
inv 2　　SN extension tableのバージョン2を表示するように指定
optyp 'amp'　　振幅を表示するように指定
bif 1;eif 1　　IF=1のデータを表示するよう指定。
inp　　パラメーターの一覧を表示して確認します。
（ここでのsnpltのパラメーターの一覧はこちら）
tvinit　　TV画面を初期化（クリア）

goと打って実行すると、TV画面上に下図のようなプロットが表示されます。縦軸の値は、SEFDの平方根です。

4.5 APCALでできたSN table をCL tableに適用する (CLCAL)

APCALによってSN extension table に格納されたSEFD（の平方根）を、CLCALによってCL extension tableに適用します。

task 'clcal'　　CLCALというtaskの使用宣言
SNVER 2　　SN extension tableのversion 2を使用
GAINVER 2　　CL extension tableの version 2を補正元に指定
gainu 3　　補正後の結果をCL extension table version 3に格納
inp　　パラメーターの一覧を表示して確認します。

（ここでのclcalのパラメーター一覧はこちら）
goと打って実行します。終了したら、imhで CL extension tableのバージョン3が追加されたことを確認しましょう。
また、imhで表示して確認もしておいてください。
以上で振幅の較正は終了です。ふう、長かったですね。では、位相の較正に行きましょー！
←前へ↑目次へ→次へ

タグ：

+ タグ編集

「A-4.ビジビリティー振幅の較正」をウィキ内検索

最終更新：2007年10月04日 20:52

添付ファイル

aips @ ウィキ