Power Spectrum을 이용한 태양 흑점 주기 계산
지난 1월 4일 게시물에서 IDL에서 Power Spectrum을 얻기 위하여 편리하게 사용할 수 있는 FFT_POWERSPECTRUM라는 함수를 소개해드렸었는데요 (링크 참조). 오늘은 이 함수의 실전 응용 사례로서, 태양의 장기간 흑점수 변화를 바탕으로 활동주기를 산출하는 작업을 해보기로 하겠습니다. 오늘 사용할 데이터는 1749년 1월부터 2015년 10월까지의 월평균 흑점수 자료로서, SIDC(Solar Influences Data Center)에서 제공되는 것을 받아서 사용하였습니다. 참고로 이 자료는 아래 링크의 웹페이지에서 Total Sunspot Number의 Monthly Mean Total Sunspot Number [1/1749 - now]에서 항상 최신 업데이트 버전을 다운로드받을 수 있습니다.
http://www.sidc.be/silso/datafiles
이 파일을 받아보면 SN_m_tot_V2.0.txt라는 이름으로 저장이 될겁니다. 그리고 그 내용을 보면 5개의 컬럼들로 구성되어 있으며, 우리의 관심사는 네번째 컬럼에 있는 월별 흑점수가 됩니다. IDL에서 이러한 파일을 읽는 가장 편리한 방법은 아무래도 IDL Astro 라이브러리에 있는 READCOL이라는 명령을 사용하는 것이라 생각되므로, 이를 사용하여 다음과 같이 읽도록 하겠습니다.
file = 'SN_m_tot_V2.0.txt'
READCOL, file, yr, mn, dfn, sn, FORMAT='I,I,F,F'
HELP, sn
여기서는 첫번째부터 네번째까지의 컬럼이 각각 year, month, fractional year, sunspot number에 해당되며, READCOL을 사용하여 각각 yr, mn, dfn, sn이라는 배열로 읽었습니다. 여기서는 sn만 사용할 생각입니다. HELP 명령을 사용하여 배열 sn에 대하여 확인해보면 3202개의 실수값들로 이루어져 있음을 알 수 있습니다. 이제 다음과 같이 이 sn 배열에 담긴 흑점수 데이터를 먼저 플롯으로 그려봅시다.
win1 = WINDOW(DIMENSIONS=[1000, 400])
pl = PLOT(sn, SYMBOL='circle', /SYM_FILLED, LINESTYLE=6, COLOR='dodger blue', $
XTITLE='Months since Jan. 1749', YTITLE='Monthly Sunspot Count', $
MARGIN=[0.08, 0.11, 0.03, 0.11], /CURRENT)
그러면 다음과 같은 그림을 얻게 됩니다.
이 그림을 그리는데 있어서는 X축을 날짜 기반으로 변환하지는 않고 그냥 1749년 1월 이후 경과한 개월수로 그렸습니다. 지금 우리가 하고자 하는 작업에 있어서는 이렇게 표출하는 것도 나쁘지 않습니다. 어쨌든 이 그림을 보면 수백년에 걸친 장기간 동안의 흑점수 변화에 있어서 주기성이 분명히 존재한다는 것이 한눈에 보입니다. 흔히 잘 알려진 약 11년의 태양활동 주기라고 하는 것은 사실상 이러한 데이터를 기반으로 얘기하는 것입니다. 그런데, 좀 더 자세히 들여다보면 11년보다 더 긴 또 다른 주기성이 존재한다는 느낌이 듭니다. 즉, 이 데이터에서는 다중의 변동주기들이 내재되어 있는 셈입니다. 원래 이와 같은 변동성 데이터를 분석하는데 있어서는자기상관(Auto-correlation)도 하나의 방법이 될 수 있겠지만, 이와 같이 다중의 주기성을 갖는 데이터를 분석하는데 있어서는 이보다는 FFT에 의한 power spectrum 분석이 더 효과적입니다. 따라서 다음과 같이 FFT_POWERSPECTRUM 함수를 사용하여 주파수 도메인에서의 그림을 그려봅시다.
f = FFT_PowerSpectrum(sn, FREQ=freq)
win2 = WINDOW(DIMENSIONS=[800, 600])
pf1 = PLOT(freq, f, /YLOG, XTITLE='Frequency (Hz)', YTITLE='Power', $
COLOR='orange red', /CURRENT, LAYOUT=[1, 2, 1])
pf2 = PLOT(freq, f, /YLOG, XTITLE='Frequency (Hz)', YTITLE='Power', $
XRANGE=[0, 0.1], COLOR='orange red', /CURRENT, LAYOUT=[1, 2, 2])
IF fig_sav THEN win.Save, 'figures/monthly_sn_power.png', WIDTH=800
특히 하단의 그림을 보면 가장 먼저 나타나는 peak의 위치는 X축 좌표가 대략 0.0076 정도로 파악이 됩니다. 이 값은 Hz 단위의 주파수입니다. 따라서 주기를 구하려면 이 값의 역수를 취해주면 됩니다. 0.0076의 역수를 구하기 위하여 1/0.0076을 계산해보면 약 131.6의 값이 나옵니다. 단 여기서 이 값의 단위는 앞서 처음에 그렸던 플롯상의 X축의 단위와 같습니다. 즉 월(month)입니다. 따라서 이를 년(year)으로 환산해보면 약 10.96이란 값이 산출됩니다. 바로 우리가 흔히 아는 태양활동주기 11년이라는 값이 이렇게 산출됩니다.
그리고 앞서 언급했듯이 또 다른 장주기 변동의 패턴이 보이는 것 같았는데, 이에 해당되는 주파수는 X축상에서 0.001 정도에 해당되는 위치로 추정이 됩니다. 이 값을 역시 월/년 단위로 환산해보면 약 83년(1000개월) 정도로 산출됩니다. 따라서 태양흑점의 장기간 갯수 변화로 살펴본 태양활동의 강도 변화는 11년 정도의 단주기성 변화외에도 약 80여년 정도의 장주기성 변화의 특성이 혼재되어 있는 모습을 보이는 것으로 해석해볼 수 있습니다. 물론 11년, 80년이 어떤 의미이고 원인이 무엇이냐의 문제는 과학적으로 풀어야 할 숙제가 되겠지요. 특히 장주기 변동 패턴의 실제 존재 여부에 대해서는 다소 의견이 분분한 측면도 있는 것도 사실입니다.
어쨌든 IDL에서 파워 스펙트럼 기법을 사용하는 방법 및 응용 사례를 살펴보았습니다. 사실 1차원 데이터외에도 이미지와 같은 2차원 데이터에 대하여 FFT 및 파워 스펙트럼 기법을 적용하려면, 원론적으로 따지면 FFT 기법 자체에 대한 이해 및 응용 방법을 알아보는 것도 의미가 있겠지만 블로그에서 다루기엔 내용이 좀 복잡한 감이 있습니다. 따라서 혹시 이러한 방법론에 관심이 있으시다면, IDL이 설치될 때 함께 제공되는 image.pdf라는 PDF 문서를 참조하시면 됩니다. 이 파일은 idl85/help/pdf 디렉토리에 있으며, 이 문서에서 7장 Transforming Between Domains라는 내용이 바로 IDL에서 FFT 기법을 사용하는 방법에 관한 내용입니다.
그리고 여기서 사용했던 데이터 파일도 아래에 첨부합니다. 2015년 10월까지의 데이터가 수록되어 있는 버전이며, 위에서 소개한 링크를 통해서 항상 최신 업데이트 버전을 받을 수 있다는 점 참조하시기 바랍니다.
SN_m_tot_V2.0.txt