太陽黑子活動(dòng)周期的分析

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

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

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

2500-1000-1500-2000-2500-3000-2000-10000RealAxis1000 2000 3000FourierCoefficientsintheComplexPlaneooo2500-1000-1500-2000-2500-3000-2000-10000RealAxis1000 2000 3000FourierCoefficientsintheComplexPlaneooo2oooooo505ooo5-功率與頻率的關(guān)系曲線則被定義為周期圖。 用周期圖法對(duì)Wolfer數(shù)進(jìn)行功率譜密度分析。程序如下：figuren=length(Y);%取FFT結(jié)果長(zhǎng)度為npower=abs(Y(1:n/2))■八2;%取前n/2個(gè)數(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為清楚起見，對(duì)功率和頻率的前50個(gè)分量作它的周期圖。程序如下: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;%將橫坐標(biāo)單位轉(zhuǎn)換為年/周期plot(period,power);%畫周期圖axis([05002e+7]);ylabel('Power');xlabel('Period(Years/Cycle)');%單位為年/周期pause

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