https://www.coursera.org/learn/data-driven-astronomy
Week 1: Thinking about data
- Principles of computational thinking
전자계산학에 기초한 자료처리의 원리
- Discovering pulsars in radio images
전파망원경 관측영상에서 펄서 찾기
----------------------------------------------
----------------------------------------------
[강의대본]
천문학에서 잡음이 깔린 자료에서 신호를 잡아내려고 시도하는 경우가 흔하다. 다음의 관측 영상에서 얼마다 많은 펄사들을 찾을 수 있을까? 머치슨 와이드필드 배열형 전파 망원경(MWA-Murchison Widefield Array)으로 관측한 영상이다.
MWA는 서호주에 위치한 저주파수 대역의 전파 망원경이다. 주파수가 80에서 300메가 헤르츠(Mhz)사이의 신호를 관측 할 수 있다. (대개 전파망원경이 수 GHz인 것에 비하면 매우 낮은 주파수 대역이다) 이 범위의 주파수는 대부분 라디오 방송국과 비슷하다. (FM방송국,아마추어 무선 통신,위성통신,항공통신 등) 이 전파 망원경의 매우 넓은 시각특성 덕분에 넓은 영역의 탐사 관측 프로젝트에 매우 유리하다. (전파 망원경의 특성:분해능은 낮지만 넓은 관측각)
이것이 MWA의 전형적인 관측 영상이다. 천체로부터 방출되는 에너지 밀도(플럭스)의 차이가 영상에서 회색의 명도 차이로 나타나 있다. 검을 수록 방출 에너지 밀도가 강하다는 듯이며 옅은 회색은 배경잡음 이다. 영상에서 보이는 검은 점들의 대부분은 멀리 떨어진 은하들로 전파를 방출 하고 있다. 몇개의 우리은하 내부의 전파방출 천체들도 보인다. 바로 펄사나 초신성 잔재들이다.
전파 천문학에서 플럭스 밀도의 측정치 단위는 잰스키(Jy)다. 면적당(미터제곱) 주파수 마다 통과하는 에너지의 량(와트)을 플럭스라 하는데 1 잰스키는 10의 마이너스 26승이다. 그러니까 플럭스 밀도란 망원경의 감지기(수신 안테나)에 입사하는 면적당 주파수별 신호의 세기(파워 스펙트럼)를 보여주는 것이다. 이에대해 좀더 알고 싶다면 이 동영상 강좌 편에 함께 제공된 자료를 읽어보기 바란다.
[참조] "마구잡이 수학", "1.4 요약 및 추가 연습문제"의 "문제 1.12 스테판-볼츠만(Stefan-Boltzmann) 법칙"
일단 이해를 돕기 위해 우리가 측정 하려는 것이 무엇인지 알려주고 싶다. 그것은 특정 주파수에서 어떤 펄사가 내는 신호의 강도(겉보기 밝기)라는 것이다. 먼저 주목할 것은 천체에서 나오는 빛(에너지 량)이 매우 미약하다는 점이다. 휴대전화에서 방출되는 전파의 세기와 비교해 보라.
천체관측 영상은 보통 FITS라는 파일의 형식으로 저장된다. 이 형식의 영상은 DS9이나 알라딘(Aladin) 같은 온-라인 도구를 통해서 볼 수 있다.
1. SAOImage DS9
2. Aladin Sky Atlas
연구 목적에 따라 관측영상을 가상의 색으로 표시하는데 전파는 색으로 표현할 수 있는 주파수가 아니라는 점을 염두에 두자. (전파는 가시광선이 아니다.) 이렇게 색깔별로 표현되어 지도의 등고선 처럼 보이는 영상은 단지 신호강도의 차이를 색으로 보여주기 위한 것이다. 지금 보여주는 영상은 우리의 커다란 전파관측 영상의 일부분을 따온 것이다.
왼쪽에 있는 영상의 중앙에 펄사가 관측된 것이 보인다. 오른쪽 영상에는 펄사가 관측되지 않았다. 우리가 뭔가 관측됐다고 말하려면, 해당 플럭스 밀도의 세기가 주변의 잡음보다 다섯배 이상의 표준편차를 가져야 한다. 만일 우리가 MWA의 관측 주파수 대역에서 전파 방출을 관측하고 있는데 펄사가 존재한다는 위치를 향하고 있음에도 대부분 관측되지 않다가 어쩌다 한번씩 관측되기도 한다.
만일 아무것도 관측하지 못했다면 매우 여러가지 원인이 있을 것이다. 펄사는 매우 먼 천체라는 것이고, 이 천체가 내뿜는 에너지 량이 엄청 나더라도 우리가 관측 하는 주파수 영역 밖에 있을 수도 있고, 깜빡이는 중간에 사진을 찍었을 수도 있다. 우리가 펄사의 신호를 감지하지 못한 정확한 이유를 대지 못할 지라도 감지하지 못했다는 사실에서 조차 뭔가 얻어낼 것이 있지 않을까?
어쨌든 천문학자들은 비록 한순간 감지하지 못할 수도 있지만 유효한 관측의 결과를 얻기 위해 현명한 방법을 개발했다. 그중 하나가 겹쳐쌓기(stacking)인데 어쩌다 감지하지 못했던 신호의 빈도를 확률적으로 찾는 방법이다.
* 이후 "겹쳐 쌓기"의 과정이 가우시안 무작위 잡음 제거(Gaussian Random Noise Reduction)효과를 설명 하는 부분에 약간 오류가 있어 보임.
[2:51] -----------------------
Stacking works because the noise in a radio image is roughly random, with a Gaussian distribution centered on zero. When you add regions of an image that just have noise, the random numbers cancel out. But when you add regions of an image in which there are signals, the signals add together, increasing what we call the signal to noise ratio.
--------------------------------
이 효과를 시각적으로 보이기 위해 1차원 신호의 예를 들어 보자.
[3:15] --------------------------
Let's take a signal that looks like a single Gaussian and then add some random noise.
단일 가우시안처럼 보이는 신호가 있다고 하자. 여기에 임의의 무작위 잡음을 더해보자.
잡음을 충분히 첨가하면 더이상 신호를 알아볼 수 없게된다.
Now say we have 100 of these signals, each with random noise added.
이제 이렇게 무작위 잡음를 첨가한 신호 100개를 준비하자.
If we take the mean of these, we can see the signal to noise ratio has increased
이 신호들을 모두 겹쳐 더한 후 평균을 구하면 신호대 잡음비가 향상된 결과를 얻을 수 있다.
and we can detect the underlying signal over the whole population of 100 sources.
100개의 원시 자료에서 빈도를 취합하여 잡음에 뭍혀있던 신호를 찾아낸 것이다.
-----------------------------------------------
이제 이 기법을 활용하여 우리가 원했던 펄사 검출 문제를 해결해 보자. 그전에 집고 넘어갈 것이 있다. 탐지되지 않은 펄사는 하늘 전역에 어디든 존재한다. 따라서 겹쳐쌓기를 하기 전에 먼저 그들의 위치를 조정해 줘야 하는데 찾고자 하는 펄사를 픽셀에 맞춰 중심에 둬야 한다. 겹쳐쌓기하는 과정은 대략 다음과 같다.
겹쳐쌓기의 평균을 계산하기 위해 영상속 전체 화소의 평균을 취한다. (한 영상 전체의 평균 밝기보정) 그결과 새로운 (정규화)영상이 만들어졌다. 다음 단계의 과정도 매 영상마다 처리해 주어야 한다. 학생들에게 주어진 과제의 끝은 여기까지다. 그리고 여기까지는 모든 과정은 아주 직관적이다.
하지만 간단해 보이는 추가사항이 하나더 주어져도 일은 아주 복잡해 질 수 있다. 추가사항 이란 전체평균(mean)대신 중간값평균(median) 겹쳐쌓기를 적용하라는 것이다. 중간값평균이 통계적으로 전체평균에 비해 매우 안정적(robust)하기 때문이다. 과학적으로 이것은 어떤 의미를 가질까? 아마 고교수학 과정에서 전체평균이 중간값평균에 비해 튀는 값에 더 쉽게 영향을 받는 다는 사실을 배웠을 것이다.
만일 대상 값의 분포가 대칭이라면 전체평균이나 중간값평균은 같게 나올 것이다.
하지만 비대칭 분포를 취하는 경우 특이값에 두 평균은 상당한 차이를 보인다. 중간값평균이 좀더 중심에 위치하게 되는 것을 볼 수 있다. 다음 강좌에서 이에 대해 좀더 자세히 살펴보기로 하겠다. 어찌 이렇게 단순해 보이는 한가지 변경으로 인해 결과를 얻기까지 계산의 수월함에 큰 영향을 주게될지 알아보기로 하자.
--------------------------------------------------
[파이썬 코딩 실습]
1. 파이썬의 배열형(Array)
#-------------------------------------------------
# Write your calculate_mean function here.
def calculate_mean(a):
return sum(a)/len(a);
# You can use this to test your function.
if __name__ == '__main__':
# Run your `calculate_mean` function with examples:
mean = calculate_mean([1, 2.2, 0.3, 3.4, 7.9])
print(mean)
mean = calculate_mean([1.2, 3.8, 2.2, 8.2, 7.1])
print(mean)
#-------------------------------------------------
#-------------------------------------------------
import numpy as np
def calc_stats(filename):
data = np.loadtxt(filename, delimiter=',')
print(data)
mean = np.mean(data)
median = np.median(data)
return np.round(mean, 1), np.round(median, 1)
# You can use this to test your function.
if __name__ == '__main__':
result = calc_stats('data.csv')
print(result)
result = calc_stats('data2.csv')
print(result)
result = calc_stats('data3.csv')
print(result)
#-------------------------------------------------
def calc_stats(filename):
data = np.loadtxt(filename, delimiter=',')
print(data)
mean = np.mean(data)
median = np.median(data)
return np.round(mean, 1), np.round(median, 1)
# You can use this to test your function.
if __name__ == '__main__':
result = calc_stats('data.csv')
print(result)
result = calc_stats('data2.csv')
print(result)
result = calc_stats('data3.csv')
print(result)
#-------------------------------------------------
#-------------------------------------------------
# Write your mean_datasets function here
import numpy as np
def mean_datasets(filename):
n = len(filename);
if (n>0):
data = np.loadtxt(filename[0], delimiter=',')
for i in range(1,n):
data += np.loadtxt(filename[i], delimiter=',')
#Mean
data_mean = data/n
return np.round(data_mean,1)
else:
print('no file')
return
# Run your function with examples from the question:
print(mean_datasets(''))
#-------------------------------------------------
import numpy as np
def mean_datasets(filename):
n = len(filename);
if (n>0):
data = np.loadtxt(filename[0], delimiter=',')
for i in range(1,n):
data += np.loadtxt(filename[i], delimiter=',')
#Mean
data_mean = data/n
return np.round(data_mean,1)
else:
print('no file')
return
# Run your function with examples from the question:
if __name__ == '__main__':
print(mean_datasets(['data1.csv', 'data2.csv', 'data3.csv']))print(mean_datasets(''))
#-------------------------------------------------
#-------------------------------------------------
# -*- coding: utf-8 -*-
"""
Created on Tue Jun 19 11:10:45 2018
@author: GoodKook
"""
# Write your load_fits function here.
from astropy.io import fits
import numpy as np
def load_fits(filename):
hdulist = fits.open(filename)
data = hdulist[0].data
arg_max = np.argmax(data)
max_pos = np.unravel_index(arg_max, data.shape)
return max_pos
if __name__ == '__main__':
szFITSname = 'image3.fits'
# Run your `load_fits` function with examples:
bright = load_fits(szFITSname)
print(bright)
# You can also confirm your result visually:
#from astropy.io import fits
import matplotlib.pyplot as plt
hdulist = fits.open(szFITSname)
data = hdulist[0].data
# Plot the 2D image data
plt.imshow(data.T, cmap=plt.cm.viridis)
plt.colorbar()
plt.show()
#-------------------------------------------------
"""
Created on Tue Jun 19 11:10:45 2018
@author: GoodKook
"""
# Write your load_fits function here.
from astropy.io import fits
import numpy as np
def load_fits(filename):
hdulist = fits.open(filename)
data = hdulist[0].data
arg_max = np.argmax(data)
max_pos = np.unravel_index(arg_max, data.shape)
return max_pos
if __name__ == '__main__':
szFITSname = 'image3.fits'
# Run your `load_fits` function with examples:
bright = load_fits(szFITSname)
print(bright)
# You can also confirm your result visually:
#from astropy.io import fits
import matplotlib.pyplot as plt
hdulist = fits.open(szFITSname)
data = hdulist[0].data
# Plot the 2D image data
plt.imshow(data.T, cmap=plt.cm.viridis)
plt.colorbar()
plt.show()
#-------------------------------------------------
4-2. 여러 FITS 형식 파일의 평균(mean) 구하기
#-------------------------------------------------
# -*- coding: utf-8 -*-
"""
Created on Tue Jun 19 12:12:13 2018
@author: GoodKook
"""
# Write your mean_fits function here:
from astropy.io import fits
import numpy as np
def mean_fits(files):
n = len(files)
if n > 0:
hdulist = fits.open(files[0])
data = hdulist[0].data
hdulist.close()
for i in range(1, n):
hdulist = fits.open(files[i])
data += hdulist[0].data
hdulist.close()
mean = data / n
return mean
if __name__ == '__main__':
# Test your function with examples from the question
data = mean_fits(['image0.fits', 'image1.fits', 'image2.fits'])
print(data[100, 100])
# You can also plot the result:
import matplotlib.pyplot as plt
plt.imshow(data.T, cmap=plt.cm.viridis)
plt.colorbar()
plt.show()
#-------------------------------------------------
"""
Created on Tue Jun 19 12:12:13 2018
@author: GoodKook
"""
# Write your mean_fits function here:
from astropy.io import fits
import numpy as np
def mean_fits(files):
n = len(files)
if n > 0:
hdulist = fits.open(files[0])
data = hdulist[0].data
hdulist.close()
for i in range(1, n):
hdulist = fits.open(files[i])
data += hdulist[0].data
hdulist.close()
mean = data / n
return mean
if __name__ == '__main__':
# Test your function with examples from the question
data = mean_fits(['image0.fits', 'image1.fits', 'image2.fits'])
print(data[100, 100])
# You can also plot the result:
import matplotlib.pyplot as plt
plt.imshow(data.T, cmap=plt.cm.viridis)
plt.colorbar()
plt.show()
#-------------------------------------------------
댓글 없음:
댓글 쓰기