太陽黑子活動周期的分析

一、引言太陽黑子是人們最早發(fā)現(xiàn)也是人們最熟悉的一種太陽表面活動。 因為太陽內(nèi)部磁場發(fā)生變化，太陽黑子的數(shù)量并不是固定的，它會隨著時間的變化而上下波動，每隔一定時間會達到一個最高點，這段時間就被稱之為一個太陽黑子周期。太陽黑子的活動呈現(xiàn)周期性變化是由施瓦貝首次發(fā)現(xiàn)的。沃爾夫 （R.Wolfer）繼而推算出11年的周期規(guī)律。實際上，太陽黑子的活動不僅呈11年的周期變化，還有海耳在研究太陽黑子磁場分布時發(fā)現(xiàn)的22年周期；格萊斯堡等人發(fā)現(xiàn)的80年周期以及蒙德極小期等。由于太陽黑子的活動規(guī)律極其復雜，時至今日科學家們?nèi)栽谂ρ芯科鋬?nèi)在的規(guī)律和特性。事實上，對太陽黑子活動規(guī)律的研究不僅具有理論意義，而且具有直接的應用需求。太陽黑子的活動呈現(xiàn)周期性變化的，沃爾夫（R.Wolfer）根據(jù)在過去的288年（1700年~1987年）間每年太陽黑子出現(xiàn)的數(shù)量和大小的觀測數(shù)據(jù)推算出11年的周期規(guī)律。我們利用Matlab強大的數(shù)據(jù)處理與仿真功能，對Wolfer數(shù)進行功率譜密度分析從而可以得到對太陽黑子活動周期的結(jié)論。實驗原理在該實驗中，對Wolfer數(shù)序列做FFT變換后得到丫（長度為n），只取其前-2個數(shù)據(jù)的功率譜密度的估計值|Y|2。原因是時域為離散的實序列的傅立葉變換對應于具有周期性且偶對稱的頻域特性，因此丫的前-個數(shù)據(jù)已經(jīng)包含了2Wolfer數(shù)的全部信息。根據(jù)DFT的頻域單位k與DTFT的頻域單位的表達關系式 k以及與f對應關系=2二丄，可以看出k與f呈線性關系N fsk;同樣地，因為丫的前-個數(shù)據(jù)已經(jīng)包含了Wolfer數(shù)的全部信息，N 2只取前-個數(shù)據(jù)分析功率-頻率圖時，對應的橫坐標時應取f二JSko（注：2N2fs為米樣頻率）

三、方案組成與流程圖1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111四、實現(xiàn)步驟及內(nèi)容讀取數(shù)據(jù)并繪制Wolfer圖，觀察太陽黑子的活動現(xiàn)象。程序如下：loadsunspot_sun.dat%讀取數(shù)據(jù)year=sunspot(:,1);%讀取年份信息wolfer=sunspot(:,2);%讀取黑子活動數(shù)據(jù)figure%新建圖像plot(year,wolfer)%畫出時域圖xlabel('Years');ylabel('SunspotData');title('SunspotData')%標注橫縱坐標和標題pauseSunspotDataYears為清楚起見，畫出部分時間段(如實驗中采取開始50年)的Wolfer圖進行觀察。

程序如下：figureplot(year(1:50),wolfer(1:50),'b.-');%用前50年數(shù)據(jù)做圖xlabel('Years');ylabel('SunspotData');title('Atthefirst50years')%標注橫縱坐標和標題140DLODSnnuAtthefirst50years1201008060140DLODSnnuAtthefirst50years1201008060402001700 1705 1710 1715 1720 1725 1730 1735 1740 1745 1750Years對已經(jīng)得到的Wolfer數(shù)應用FFT技術分析它的變化規(guī)律。程序如下：Y=fft(wolfer);%對全部數(shù)據(jù)做FFTY(1)=[];%由于算法問題，Y⑴為所有fft數(shù)值之和，因此需舍棄第一個點pause觀察Wolfer數(shù)在FFT分析后得到的在復平面中的分布圖。程序如下：figureplot(Y,'ro')%在復平面做圖，空心點title('FourierCoefficientsintheComplexPlane');xlabel('RealAxis');%標注橫坐標為’實部’ylabel('ImaginaryAxis');%標注縱坐標為’虛部’pause

2500-1000-1500-2000-2500-3000-2000-10000RealAxis1000 2000 3000FourierCoefficientsintheComplexPlaneooo2500-1000-1500-2000-2500-3000-2000-10000RealAxis1000 2000 3000FourierCoefficientsintheComplexPlaneooo2oooooo505ooo5-功率與頻率的關系曲線則被定義為周期圖。 用周期圖法對Wolfer數(shù)進行功率譜密度分析。程序如下：figuren=length(Y);%取FFT結(jié)果長度為npower=abs(Y(1:n/2))■八2;%取前n/2個數(shù)據(jù)求其功率nyquist=1/2;%取最大頻率為0.5freq=(1:n/2)/(n/2)*nyquist;%將FFT轉(zhuǎn)換為DTFTplot(freq,power)%畫周期圖xlabel('cycles/year');title('Periodogram')pause77為清楚起見，對功率和頻率的前50個分量作它的周期圖。程序如下:figureplot(freq(1:50),power(1:50))%為觀察方便，取前50年數(shù)據(jù)做圖xlabel('cycles/year')pause

21.81.61.41.210.80.60.40.20X100 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18cycles/yearX100 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18cycles/yearfigureperiod=1./freq;%將橫坐標單位轉(zhuǎn)換為年/周期plot(period,power);%畫周期圖axis([05002e+7]);ylabel('Power');xlabel('Period(Years/Cycle)');%單位為年/周期pause

(5)在功率與周期關系曲線圖中標出功率的最高點，該位置對應的周期則為太陽黑子活動周期。holdon;index=find(power==max(power));%找到頻率最大點，該點橫坐標即為太陽黑子周期mainPeriodStr=num2str(period(index));plot(period(index),power(index),'r.','MarkerSize',25);%用實心點指出該點text(period(index)+2,power(index),['Period=；mainPeriodStr]);%文字標注該點holdoff;7272X101.81.61.41.2w10.80.60.40.20?Period=?Period=11.0385Period(Years/