데이터분석 — 상관처리 데이터로 영상 만들기

관측데이터의 상관처리에서 살펴본 상관기의 출력은 아직 우리가 보는 "사진"이 아니라, 각 망원경 쌍(기선)·시간·주파수별로 얻은 밝기와 위상 정보인 가시성 데이터(visibility data)입니다. 상관처리가 "신호를 맞추어 합치는" 과정이라면, 이 가시성 데이터를 실제로 눈에 보이는 영상으로 바꾸는 것은 완전히 다른 별도의 계산 과정이며, 이 페이지에서 그 방법을 살펴봅니다. 자료는 이 과정에서 상관기가 만든 FITS-IDI → 보정을 마친 UVFITS → 완성된 FITS 영상으로 모습을 바꿉니다. 먼저 이 자료들을 담는 FITS 파일의 구조를 보고, 처리 단계와 영상이 만들어지는 원리(역푸리에 변환과 CLEAN), 화면에 표시하는 방법, 파이썬 분석과 연구자의 고민까지 차례로 따라가 봅니다.

FITS 파일 — 천문 자료를 담는 표준 그릇

FITS는 1979년에 천문대끼리 자료를 주고받기 위해 만든 형식으로, 지금은 전파·광학·X선 등 거의 모든 천문 자료가 이 형식을 씁니다. 구조는 단순합니다. 파일은 2,880바이트 블록을 이어 붙인 것이고, 그 안에 헤더(설명서)와 데이터(숫자 배열)가 한 쌍을 이룬 HDU(Header-Data Unit)가 하나 이상 들어 있습니다. 헤더는 사람이 그대로 읽을 수 있는 80글자짜리 줄(카드)들이고, 데이터는 컴퓨터가 읽는 이진 숫자입니다.

전파 영상 FITS에서 특히 중요한 헤더 키워드는 다음과 같습니다. 헤더만 읽어도 "어디를, 무엇으로, 어떤 크기로 찍었는지" 알 수 있습니다.

VLBI 자료의 세 가지 FITS

같은 FITS라도 처리 단계마다 담긴 내용이 다릅니다. 앞의 두 가지는 아직 '그림'이 아니라 측정값의 표이고, 마지막 하나가 우리가 보는 영상입니다.

파일 안을 열어 보면 — FITS-IDI와 UVFITS의 구조

FITS-IDI는 내용이 빈 첫 HDU 뒤에 이진 테이블 HDU가 여러 개 이어진 파일입니다. 엑셀의 시트 여러 장과 비슷하게, 시트마다 이름(EXTNAME)이 붙어 있습니다.

UVFITS는 오래된 '랜덤 그룹' 구조를 씁니다. 측정 하나가 그룹 하나이고, 그룹마다 앞에 매개변수(UU, VV, WW, BASELINE, DATE 등)가 붙고 뒤에 (실수, 허수, 가중치) 묶음이 편파·채널·대역별로 이어집니다. 헤더의 PTYPEn이 매개변수 이름을, CTYPEn이 COMPLEX·STOKES·FREQ·IF·RA·DEC 축 이름을 알려 줍니다. 안테나 목록과 주파수 정보는 'AIPS AN', 'AIPS FQ' 같은 이름의 테이블 HDU로 뒤에 붙습니다.

FITS-IDI에서 FITS 영상까지 — 처리 단계 따라가기

상관기 출력을 영상으로 만드는 일은 크게 보정 → 역푸리에 변환 → CLEAN입니다. 주로 미국 국립전파천문대(NRAO)가 만든 AIPS와 최신 표준 소프트웨어 CASA, VLBI 영상화에 특화된 Difmap을 쓰며, KVN·EAVN 자료도 이런 프로그램으로 처리합니다. 실제 작업 단위로 나누면 다음과 같고, 이름은 달라도 하는 일은 같습니다.

7단계 자세히 — 역푸리에 변환과 CLEAN

가시성 데이터는 하늘의 밝기 분포를 푸리에 변환한 값(공간주파수)을 u-v 평면의 여러 점에서 잰 것입니다. 기선 하나는 한 순간에 u-v 평면의 점 하나(와 그 맞은편 점)를 재고, 지구가 자전하면 기선의 길이와 방향이 바뀌어 점들이 타원 모양의 궤적을 그립니다(간섭계 원리). 이 점들을 격자에 올려 역푸리에 변환하면 영상이 되지만, 망원경 수가 유한해 u-v 평면에 빈칸이 많으므로 진짜 천체에 가짜 무늬(사이드로브)가 겹친 보정 전 영상(dirty image)이 나옵니다. 점 하나를 관측했을 때 생기는 이 무늬를 보정 전 빔(dirty beam)이라고 합니다.

CLEAN은 이 무늬를 걷어 내는 대표적인 디컨볼루션 방법입니다. ① 영상에서 가장 밝은 곳을 찾아 ② 그 밝기의 일부(보통 10 %)만큼 보정 전 빔 무늬를 빼고 ③ 그 위치와 세기를 'CLEAN 성분'으로 기록하는 일을 수백 ~ 수천 번 되풀이합니다. 마지막으로 기록한 성분들을 사이드로브가 없는 깨끗한 타원(복원 빔, FITS 헤더의 BMAJ·BMIN·BPA)으로 다시 그리고, 빼고 남은 잔차를 더하면 최종 영상이 됩니다. 이것이 FITS 영상으로 저장됩니다.

👉 망원경 구성과 관측 시간을 바꿔 u-v 채움과 보정 전 빔이 어떻게 달라지는지 보고, CLEAN 버튼으로 가짜 무늬를 걷어 내 보세요!

KVN 4대로 적도 근처의 3C 273을 보면 u-v 궤적이 거의 가로줄이 되어, 보정 전 빔에 진짜 봉우리만큼 높은 가짜 봉우리가 위아래로 되풀이됩니다. 이럴 때는 CLEAN도 성분을 엉뚱한 곳에 놓기 쉽습니다(④의 초록 + 표시와 비교해 보세요). 적위가 높은 3C 84를 고르면 궤적이 타원을 그리며 u-v 평면을 더 고르게 채워 무늬가 줄어듭니다. 관측 시간을 늘려도 궤적이 길어져 무늬가 줄어듭니다. 일본 VERA를 더하면 기선이 길어져 빔(해상도)이 작아지고 u-v 채움도 촘촘해집니다. 이것이 한국·일본·중국 망원경을 묶는 동아시아 VLBI 관측망(EAVN)의 장점입니다. 처음 몇 번의 CLEAN에서 영상이 크게 바뀌고, 뒤로 갈수록 바뀌는 양이 줄어드는 것도 확인해 보세요. 밝은 성분부터 차례로 빠지기 때문입니다.

VLBI 영상화에는 Difmap을 많이 씁니다. 6단계에서 만든 UVFITS로 영상을 만드는 짧은 예는 다음과 같습니다(값은 43 GHz 관측을 가정한 예시이며, 실제로는 자료를 보며 조정합니다).

! Difmap 예시 — UVFITS → FITS 영상
observe 3C273_Q.uvf      ! 보정된 가시성 불러오기
select i                 ! 전체 세기(Stokes I)만
mapsize 256, 0.1         ! 256 × 256 픽셀, 한 칸 0.1 mas
uvweight 2, -1           ! 가중치 방식(자연·균일 사이)
clean 200, 0.05          ! CLEAN 200번, 한 번에 5 %씩
selfcal                  ! 모형으로 위상 자기보정
clean 500, 0.05
selfcal true, true       ! 진폭까지 자기보정
restore                  ! 복원 빔으로 다시 그리기
wmap 3C273_Q.fits        ! FITS 영상으로 저장

# CASA 예시 (파이썬) — 보정된 자료 → FITS 영상
tclean(vis='kv26a_3c273.ms',
       imagename='3c273_q',
       imsize=256,
       cell='0.1mas',          # 한 칸 크기
       weighting='briggs', robust=0,
       niter=2000,             # CLEAN 최대 반복
       threshold='3mJy')       # 이 세기까지 CLEAN
exportfits(imagename='3c273_q.image',
           fitsimage='3c273_q.fits')

FITS 영상 파일을 화면의 그림으로 — 표시 과정

FITS 영상 파일에는 '색'이 없습니다. 들어 있는 것은 픽셀마다의 세기(Jy/beam) 숫자뿐이고, 이것을 사람이 볼 수 있는 그림으로 바꾸는 것은 보는 프로그램의 몫입니다. 같은 FITS라도 표시 방법에 따라 전혀 다른 느낌의 그림이 됩니다.

이렇게 표시한 실제 KVN·KaVA 영상에서 축·빔·등고선·색·크기 막대를 찾아 읽어 보는 연습은 관측영상 편에서 이어집니다.

🔬 FITS 실험실 — 파일을 열어 숫자에서 그림까지

👉 모의 FITS 파일의 헤더와 바이트를 살펴보고, 스트레치·색 지도·등고선을 바꿔 그림이 어떻게 달라지는지 보세요. 내 FITS 파일을 열어 볼 수도 있습니다!

파이썬으로 FITS 다루기

영상을 만든 뒤의 분석(세기 재기, 그림 그리기, 여러 관측 비교)은 대부분 파이썬으로 합니다. 천문학 공동체가 함께 만드는 Astropy가 FITS 읽기·쓰기와 하늘 좌표 계산의 표준이고, CASA도 파이썬 모듈로 불러 쓸 수 있습니다.

예제 1 — FITS 영상 읽고 좌표·등고선·빔이 있는 그림 그리기

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from astropy.io import fits
from astropy.wcs import WCS
from astropy.visualization import simple_norm

with fits.open('3C273_Q.fits') as hdul:
    hdul.info()                          # HDU 목록 보기
    hdr = hdul[0].header
    img = hdul[0].data.squeeze()         # (1, 1, 256, 256) → (256, 256): 주파수·편파 축 제거

print(hdr['OBJECT'], hdr['BUNIT'], hdr['CRVAL3'] / 1e9, 'GHz')
pix_mas = abs(hdr['CDELT1']) * 3.6e6     # 도 → 밀리각초
bmaj, bmin = hdr['BMAJ'] * 3.6e6, hdr['BMIN'] * 3.6e6

rms = np.nanstd(img[:60, :60])           # 천체가 없는 구석으로 잡음 추정
peak = np.nanmax(img)
levels = 3 * rms * 2.0 ** np.arange(12)  # 3σ부터 2배씩 커지는 등고선
levels = levels[levels < peak]

wcs = WCS(hdr).celestial                 # 하늘 좌표 두 축만 사용
ax = plt.subplot(projection=wcs)
ax.imshow(img, origin='lower', cmap='inferno',          # FITS는 아랫줄부터 → origin='lower'
          norm=simple_norm(img, 'log', vmin=rms, vmax=peak))   # astropy 6.1 이상
ax.contour(img, levels=levels, colors='white', linewidths=0.6)
ax.add_patch(Ellipse((20, 20), bmin / pix_mas, bmaj / pix_mas, angle=hdr['BPA'],
                     facecolor='none', edgecolor='cyan'))   # 왼쪽 아래 빔 타원
ax.set_xlabel('적경'); ax.set_ylabel('적위')
plt.savefig('3C273_Q.png', dpi=200)

예제 2 — 과학 값 계산하기 (잡음, 다이내믹 레인지, 전체 세기)

beam_pix = np.pi * hdr['BMAJ'] * hdr['BMIN'] / (4 * np.log(2)) / (abs(hdr['CDELT1']) * abs(hdr['CDELT2']))
mask = img > 3 * rms                     # 3σ보다 밝은 곳만 천체로 봄
total_jy = img[mask].sum() / beam_pix    # Jy/beam의 합 ÷ 빔 넓이(픽셀) = Jy
print(f'잡음 {rms*1e3:.2f} mJy/beam, 최대 {peak:.3f} Jy/beam, '
      f'다이내믹 레인지 {peak/rms:.0f}, 전체 세기 {total_jy:.3f} Jy')

y, x = np.unravel_index(np.nanargmax(img), img.shape)
sky = wcs.pixel_to_world(x, y)           # 가장 밝은 픽셀의 하늘 좌표
print(sky.to_string('hmsdms'))

예제 3 — 영상 전 단계의 자료 열어 보기 (UVFITS, FITS-IDI)

with fits.open('3C273_Q.uvf') as uvf:     # UVFITS: 랜덤 그룹
    g = uvf[0].data
    print(g.parnames)                     # ['UU', 'VV', 'WW', 'BASELINE', 'DATE', ...] (프로그램마다 이름이 조금 다름)
    freq = uvf[0].header['CRVAL4']        # 기준 주파수(Hz) — 어느 축인지는 CTYPEn으로 확인
    u_lambda = g.par(0) * freq            # 초 × Hz = 파장 단위의 u
    vis = g.data                          # (측정 수, …, 편파, 3) 마지막 축 = 실수·허수·가중치

with fits.open('kv26a.idifits') as idi:   # FITS-IDI: 이진 테이블 여러 개
    idi.info()                            # ARRAY_GEOMETRY, ANTENNA, FREQUENCY, SOURCE, UV_DATA …
    uv = idi['UV_DATA'].data
    a1, a2 = divmod(uv['BASELINE'], 256)  # 기선 번호 → 안테나 두 개
    print(idi['ANTENNA'].data['ANNAME'])  # 안테나 이름

파일 이름은 예시입니다. 실제 키워드와 축 순서는 만든 프로그램마다 조금씩 다르므로, 먼저 hdul.info()와 헤더를 출력해 확인하는 습관이 중요합니다.

연구자는 무엇을 고민할까?

버튼 몇 번으로 영상이 나오는 것처럼 보이지만, 같은 자료로도 처리 선택에 따라 다른 영상이 나옵니다. 연구자는 단계마다 다음과 같은 질문을 붙들고 씨름하며, "이 구조가 정말 천체에 있는가, 처리 과정이 만든 것인가"를 끝까지 확인합니다.

본 교육자료는 Anthropic(Claude) 및 OpenAI(ChatGPT)를 활용하여 제작하였으며, 제작자가 내용을 검토·편집하였습니다. 주요 참고자료와 출처는 각 페이지의 「참고 자료 및 출처」에 명시하였습니다. 출처가 별도로 표기되지 않은 페이지와 참고자료에 포함되지 않은 내용은 제작자가 직접 작성하거나 정리한 것입니다. 내용에 오류가 있거나 수정이 필요한 사항을 발견하시면 kvnsite@kasi.re.kr로 알려 주시기 바랍니다.