地球物理學(xué)畢業(yè)論文

題 目 基于臨界平衡假設(shè)的 太陽(yáng)風(fēng)湍流各向異性的觀測(cè)分析 英文題目 Observation the solar wind turbulence anisotropy based on the critical balance assumption 摘 要 Elsasser 于 1950 年對(duì)磁流體方程引入新變量后 Goldreich 和 Shirida 于 1995年提出臨界平衡理論假設(shè)。并將理論應(yīng)用于太陽(yáng)風(fēng)湍流中得到功率譜各向異性特征和功率譜的譜指標(biāo)特征。 本文的主要目的是做出速度場(chǎng)的功率譜，進(jìn)一步檢驗(yàn)臨界平衡理論假設(shè)的正確性與適用性。文中采用了安裝在 WIND 衛(wèi)星上的磁場(chǎng)探測(cè)器 MFI 和三維等離子體探測(cè)器 3DP 在 1995 年 1 月 31 日所采集到的總磁場(chǎng)和等離子體速度場(chǎng)數(shù)據(jù)。并且采用基于小波變換的功 率譜分析，做出速度場(chǎng)功率譜密度函數(shù)圖和磁場(chǎng)功率譜密度函數(shù)圖。其中橫軸為局部磁場(chǎng)和太陽(yáng)徑向向外方向的夾角；縱軸為時(shí)間尺度。通過(guò)功率譜密度函數(shù)圖的分析可以得出功率譜隨角度變化的分布特征。 通過(guò)計(jì)算機(jī)的計(jì)算結(jié)果表明：當(dāng)天的速度場(chǎng)功率譜顯式出各向同性的特征，同時(shí)磁場(chǎng)的功率譜也同樣顯示出了各向同性的特征。這說(shuō)明臨界平衡理論假設(shè)并不是普遍適用的。 從原始數(shù)據(jù)看，由于速度擾動(dòng)和速度化磁場(chǎng)擾動(dòng)高度相關(guān)，筆者認(rèn)為這可能是產(chǎn)生各向同性的原因之一。遺憾的是本文只針對(duì) 1995 年 1 月 31 日的數(shù)據(jù)進(jìn)行了考察，無(wú)法同其他明顯存在阿爾文波 時(shí)間段的功率譜進(jìn)行比較。同時(shí)當(dāng)功率譜為各向同性時(shí)，功率譜譜指標(biāo)特征也有待進(jìn)一步探究。 關(guān)鍵詞 ： 臨界平衡；太陽(yáng)風(fēng)湍流；功率譜；小波分析；各向異性； Abstract Goldreich and Shirida in 1995, propose a critical balance of theoretical assumptions. After Elsasser interoduce a new varibles in the MHD equation in 1950. The theory can be applied to the solar wind turbulence characteristics of the anisotropy power spectrum and characteristics of the power spectrum spectral index. The main purpose of this paper is to make the power spectrum of the velocity field, Further test the accuracy and applicability of the critical balance of theoretical assumptions. This paper uses a magnetic field and plasma velocity field data on January 31, 1995. collected by MFI and 3DP which installed on the WIND spacecraft. And based on wavelet transform power spectrum analysis, to make the power spectral density of the velocity field maps and magnetic field power spectral density map. Where the horizontal axis is the angle between the local magnetic field and solar radial outward direction, the longitudinal axis is the time scale. By the power spectral density maps can be easily observed power spectrum changes with the angle distribution characteristics. Computational results show that day the velocity field power spectrum shows isotropic characteristics of the magnetic field power spectrum also shows isotropic characteristics. This shows that the theoretical assumptions of the critical balance is not generally applicable. From the data, highly correlated to the velocity disturbance and the speed of the magnetic field disturbance, I think this may be one of the reasons to produce isotropic. Unfortunately, this paper only examined data for January 31, 1995, can not be compared with other power spectrum when the Alfven waves is apparent existence. And when the power spectrum is isotropic, the characteristics of the power spectrum spectral index also needs further exploration. Key words: Critical balance； solar wind turbulence； power spectrum； wavelet analysis； anisotropy； 目 錄 前 言 . 錯(cuò)誤 !未定義書(shū)簽。 1 背景介紹 . 錯(cuò)誤 !未定義書(shū)簽。 1 1 太陽(yáng)大氣與太陽(yáng)風(fēng) . 錯(cuò)誤 !未定義書(shū)簽。 1 2 太陽(yáng)風(fēng)及湍流 . 3 1 3 阿爾文波 . 4 1 4 臨界平衡 . 5 1 5 太陽(yáng)風(fēng)湍流功率譜的各向異性 . 8 1 6 衛(wèi)星介紹 . . 10 2 數(shù)據(jù)及數(shù)據(jù)處理原理 . 1 錯(cuò)誤 !未定義書(shū)簽。 2 1 數(shù)據(jù)來(lái)源及處理過(guò)程 . 13 2 2 小波分析數(shù)據(jù)原理 . 13 2 2 1 小波分析原理 . 13 2 2 2 實(shí)際工作中小波分析的運(yùn)用 . 14 2 2 3 功率譜密度 PSD . 15 2 3 功率譜隨角度的分布 . 16 3 數(shù)據(jù)處理及處理結(jié)果 . 錯(cuò)誤 !未定義書(shū)簽。 8 3 1 太陽(yáng)風(fēng)中的時(shí)間序列 . 錯(cuò)誤 !未定義書(shū)簽。 8 3 2 磁場(chǎng)數(shù)據(jù)處理 . 錯(cuò)誤 !未定義書(shū)簽。 9 3 3 速度場(chǎng)數(shù)據(jù)處理 . 20 結(jié) 論 . 21 致 謝 . 22 參考文獻(xiàn) . 23 附 錄 一 . 24 附 錄 二 . 32 附 錄 三 . 36 1 前 言 1950年， Elsasser對(duì)磁流體方程進(jìn)行了改進(jìn)，引入了 Elsasser新變量。新變量將磁流體的磁場(chǎng)和動(dòng)量結(jié)合在了一起，不僅將磁流體的動(dòng)量方程和磁 感應(yīng)方程結(jié)合在一起，同時(shí)新方程“顯式”的描述了反向傳播的阿爾文波的非線(xiàn)性相互作用。 在此基礎(chǔ)上 Goldreich 和 Shirida 于 1995年提出了臨界平衡理論。他們假設(shè)磁流體中波動(dòng)的傳播效應(yīng)和反向波動(dòng)相互作用效應(yīng)相當(dāng)。并且將理論應(yīng)用于太陽(yáng)風(fēng)中，得到了太陽(yáng)風(fēng)湍流各向異性的結(jié)論。理論也給出了垂直于磁場(chǎng)方向和平行于磁場(chǎng)方向的功率譜的譜指數(shù)分別是 -5/3和 -2 。 Timothy S. Horbury于 2008年的一篇文章中給出了垂直于磁場(chǎng)方向和平行于磁場(chǎng)方向的功率譜的譜指數(shù)的觀測(cè)值，結(jié)果與理論符合的很好。 J. J. Podesta于 2009年的一篇文章中進(jìn)一步細(xì)致的研究了局部磁場(chǎng)方向與太陽(yáng)徑向向外方向不同夾角不同時(shí)的功率譜譜指標(biāo)特征，結(jié)果同樣與理論符合的很好。在 He&Marsch&Tu&Yao&Tian等人共同發(fā)表的一篇文章中，通過(guò)Helios、 STEREO飛船的觀測(cè)結(jié)果，發(fā)現(xiàn)太陽(yáng)風(fēng)湍流功率譜的各向異性特征。 很明顯，現(xiàn)有的觀測(cè)結(jié)果普遍支持臨界平衡理論。不僅觀測(cè)到了太陽(yáng)風(fēng)功率譜的各向異性特征。同時(shí)，功率譜譜指數(shù)與角度的關(guān)系也與理論吻合的很好。但是，臨界平衡這個(gè)理論現(xiàn)在還處在一個(gè)假設(shè)性的階段，并未完全得到肯定。 而且關(guān)于速度場(chǎng)的功率譜的研究工作比較少。所以本文的目的在于做出速度場(chǎng)的功率譜，進(jìn)一步檢驗(yàn)臨界平衡理論假設(shè)的正確性與適用性。這一點(diǎn)對(duì)于一個(gè)理論的發(fā)展是有幫助的。無(wú)論最后得到肯定的還是否定的結(jié)論，都不僅幫助我們改善我們現(xiàn)有的理論，而且將進(jìn)一步加深我們對(duì)于太陽(yáng)風(fēng)湍流特征的認(rèn)識(shí)。從而更好的理解日地之間的空間物理特性。 2 1 背景介紹 1.1 太陽(yáng)及太陽(yáng)大氣 太陽(yáng) 1是太陽(yáng)系的主導(dǎo)恒星，是維持人類(lèi)生存活動(dòng)的主要源泉。其質(zhì)量大約為302 10 kg，半徑是 70萬(wàn)千米。日地直線(xiàn)距離是 1.5億千米，相當(dāng)于 215個(gè)太陽(yáng)半徑。太陽(yáng)等離子體主要由氫 (約 90%)、氦 (約 10%)以及碳、氮、氧等其他元素 (約 0.1%)組成 章振大 , 1992。太陽(yáng)的電磁和粒子輻射是影響日地空間環(huán)境的主要因素。太陽(yáng)的電磁輻射覆蓋了電磁波譜中從射線(xiàn)到射電波段的極寬的頻率范圍。同時(shí)，太陽(yáng)連續(xù)不斷地向外發(fā)出數(shù)十萬(wàn)度高溫的稀薄的磁化等離子體，即以每秒數(shù)百公里的速度向外運(yùn)行的太陽(yáng)風(fēng)。 太陽(yáng)內(nèi)部從里往外大致可以分為日核、輻射區(qū)、對(duì)流層。日核 (日心到 0.25Rs)是產(chǎn)能區(qū)，太陽(yáng)在自身引力作用下，物質(zhì)向日核聚集，導(dǎo)致極高的溫度 (107K)，從而不停地進(jìn)行著由氫核聚變成氦核的熱核反應(yīng)。日核產(chǎn)生的能量在輻射區(qū) (0.25-0.75 Rs)經(jīng)過(guò)多次吸收、散射和再發(fā)射逐漸向外傳輸。輻射區(qū)之上是對(duì)流層 (0.75 Rs到太陽(yáng)表面 )，其中的物質(zhì)處于劇烈的對(duì)流狀態(tài)。 對(duì)流層之上便是我們?nèi)庋劭梢?jiàn)的太陽(yáng)表面，稱(chēng)為光球。光球及其上的色球、過(guò)渡區(qū)和日冕合稱(chēng)為太陽(yáng)大氣。圖 1.1 中的兩條曲線(xiàn)是根據(jù)著名的太陽(yáng)大氣 VAL-3C模型 Vernazza et al., 1981計(jì)算得到的太陽(yáng)大氣中溫度和粒子數(shù)密度隨高度變化的曲線(xiàn)。通常將太陽(yáng)大氣中溫度極小處作為光球和色 球的分界點(diǎn)。占太陽(yáng)輻射能量中絕大部分的可見(jiàn)光主要來(lái)源于光球。對(duì)流層的對(duì)流運(yùn)動(dòng)反映在光球上，主要表現(xiàn)為米粒和超米粒的形式。 光球上方是色球，其溫度從 4200 K上升到 20000 K左右。色球輻射主要呈現(xiàn)出亮的或暗的網(wǎng)絡(luò)結(jié)構(gòu)，網(wǎng)絡(luò)之間的區(qū)域稱(chēng)為網(wǎng)絡(luò)內(nèi)區(qū)。在高色球中存在著許多針狀的結(jié)構(gòu)，被認(rèn)為是冷而密的色球物質(zhì)沿磁力線(xiàn)向日冕噴射所形成的。 過(guò)渡區(qū)是色球之上溫度大致在 104 K和 106 K之間的區(qū)域。由圖 1-1可以看出，在色球和日冕之間的很窄的過(guò)渡區(qū)內(nèi)，太陽(yáng)大氣的溫度陡增，密度陡降。經(jīng)過(guò)數(shù)十年的研究，人們已經(jīng)認(rèn)識(shí)到 過(guò)渡區(qū)遠(yuǎn)非一個(gè)靜態(tài)的分層結(jié)構(gòu)，而是一個(gè)磁場(chǎng)和等離子體結(jié)構(gòu)非常不均勻的動(dòng)態(tài)區(qū)域。 日冕是太陽(yáng)大氣最外面的一層稀薄的等離子體。 610 K以上的高溫導(dǎo)致日冕氣體高度電離。日冕輻射主要由如下幾部分組成：自由電子對(duì)光球輻射的湯姆遜散射形成的 K冕，擴(kuò)展日冕 (2.5 Rs以上 )中的塵埃散射光球輻射而形成的 F冕，由離子譜線(xiàn)輻射組成的 E冕，以及主要由軔致輻射所產(chǎn)生的 X冕 (軟 X射線(xiàn) )和 R冕 (射電 )。日冕中的結(jié)構(gòu)有輻射很弱的冕洞、活動(dòng)區(qū)里的冕環(huán)、極區(qū)外呈羽毛狀的極羽、以及赤道附近的盔狀冕流等。 3 圖 1-1 太陽(yáng) 大氣中的溫度和粒子數(shù)密度隨高度的平均變化特征。引自 Peter 2004。 1.2 太陽(yáng)風(fēng)及湍流 太陽(yáng)風(fēng)是從太陽(yáng)上層大氣射出的超高速等離子體（帶電粒子）流。在不是太陽(yáng)的情況下，這種帶電粒子流也常稱(chēng)為“恒星風(fēng)”。 在太陽(yáng)的日冕層的高溫（幾百萬(wàn)開(kāi)氏度）下，氫、氦等原子已經(jīng)被電離成帶正電的質(zhì)子、氦原子核和帶負(fù)電的自由電子等。這些帶電粒子運(yùn)動(dòng)速度極快，以致不斷有帶電的粒子掙脫太陽(yáng)的引力束縛，射向太陽(yáng)的外圍，形成太陽(yáng)風(fēng)。 太陽(yáng)風(fēng)的速度一般在 200-800km/s。 一般認(rèn)為在太陽(yáng)極小期，從太陽(yáng)的磁場(chǎng)極地附近吹出 的是高速太陽(yáng)風(fēng)，從太陽(yáng)的磁場(chǎng)赤道附近吹出的是低速太陽(yáng)風(fēng)。太陽(yáng)的磁場(chǎng)的活動(dòng)性是會(huì)變化的，周期大約為 11年。 高溫導(dǎo)致日冕氣體膨脹，連續(xù)不斷地向外發(fā)射等離子體流，到達(dá)數(shù)個(gè)太陽(yáng)半徑的距離后變成超聲速的流動(dòng)，形成太陽(yáng)風(fēng)。太陽(yáng)風(fēng)主要由電子和質(zhì)子組成，另有少量粒子和極少量重離子。太陽(yáng)風(fēng)把太陽(yáng)磁場(chǎng)帶到行星際空間，形成螺旋狀的行星際磁場(chǎng)。 Parker 1958預(yù)言太陽(yáng)風(fēng)后不久，前蘇聯(lián)和美國(guó)的空間飛船探測(cè)即證實(shí)了太陽(yáng)風(fēng)的存在。根據(jù)地球軌道 (1 AU)和行星際空間實(shí)地測(cè)量的流速，可將太陽(yáng)風(fēng)分為高速流和低速流兩種。 太陽(yáng)風(fēng) 一詞是在 1950年代被尤金派克提出。但是直到 1960年代才證實(shí)了它的存在。長(zhǎng)期觀測(cè)發(fā)現(xiàn)，當(dāng)太陽(yáng)存在冕洞時(shí)，地球附近就能觀測(cè)到高速的太陽(yáng)風(fēng)。因此天文學(xué)家認(rèn)為高速太陽(yáng)風(fēng)的產(chǎn)生與冕洞有密切的關(guān)系。太陽(yáng)表面的磁場(chǎng)及等離子體活動(dòng)對(duì)地球有很重要的影響。當(dāng)太陽(yáng)發(fā)生強(qiáng)烈的活動(dòng)時(shí)，大量的帶電粒子隨著太陽(yáng)風(fēng)吹向地球的兩極，就會(huì)在兩極的電離層引發(fā)美麗的極光。 表 1-1列舉了 1 AU處觀測(cè)到的高速流和低速流的平均特性。 4 表 1-1 1 AU處觀測(cè)到的太陽(yáng)風(fēng)高速流和低速流的平均特性。 引自 Axford and McKenzie1997， Xia 2003和 Holzer 2005。 湍流 2，也稱(chēng)為紊流，是流體的一種流動(dòng)狀態(tài)。當(dāng)流速很小時(shí)，流體分層流動(dòng)，互不混合，稱(chēng)為層流，或稱(chēng)為片流；逐漸增加流速，流體的流線(xiàn)開(kāi)始出現(xiàn)波浪狀的擺動(dòng)，擺動(dòng)的頻率及振幅隨流速的增加而增加，此種流況稱(chēng)為過(guò)渡流；當(dāng)流速增加到很大時(shí)，流線(xiàn)不再清楚可辨，流場(chǎng)中有許多小漩渦，稱(chēng)為湍流，又稱(chēng)為亂流、擾流或紊流。 而湍流在太陽(yáng)上和日球?qū)邮呛芷毡榈默F(xiàn)象。太陽(yáng)風(fēng)也是動(dòng)蕩的，并不是平靜的流動(dòng)?？梢岳酶道锶~譜分析或是小波譜分析來(lái)研究太陽(yáng)風(fēng)中的湍流 3，來(lái)計(jì)算其磁場(chǎng)的功率譜 4，結(jié)果與早已充分發(fā)展的磁流體湍流的結(jié)果類(lèi)似。特別是，功率譜的譜指數(shù)接近 -5/35。 1.3 阿爾文波 阿爾文波 6，又稱(chēng)剪切阿爾文波，是等離子體中的一種沿磁場(chǎng)方向傳播的波，這種波的頻率遠(yuǎn)低于等離子體的回旋頻率，是一種線(xiàn)偏振的低頻橫波。處在磁場(chǎng)中的導(dǎo)電流體在垂直于磁場(chǎng)的方向上受到局部擾動(dòng)時(shí)，沿著磁感線(xiàn)方向的磁張力提供恢復(fù)力，就會(huì)激發(fā)阿爾文波。阿爾文波是由瑞典物理學(xué)家漢尼斯阿爾文首先預(yù)言的，因此得名。后來(lái)隆德奎斯特使用 1特斯拉左右的磁場(chǎng)在水銀中觀察到了阿爾文波，列納爾特使用 液態(tài)鈉也證實(shí)了阿爾文波的存在。阿爾文波的色散關(guān)系為： 這樣的波稱(chēng)為斜阿爾文波。 =0 時(shí)是沿著磁感線(xiàn)的方向傳播的，稱(chēng)為阿爾文波。此時(shí)阿爾文波的相速度為： 稱(chēng)為斜傳阿爾文速度，其中 0 是等離子體的磁導(dǎo)率， 0是離子的密度。在垂直于磁感線(xiàn)的方向上阿爾文波不能傳播。 5 1.4 臨界平衡 關(guān)于太陽(yáng)風(fēng)的傳播， 1995 年 Goldreich 和 Shirida 提出基于臨界平衡假設(shè)的理論，這里做一個(gè)簡(jiǎn)單的介紹。 我們先來(lái)看一下非線(xiàn)性簡(jiǎn)化磁流體 (MHD)方程的新形式 7。 首先，動(dòng)量方程是： 2( ) ( )t o tu u u P v u b bt (1.1) 這個(gè)方程將速度進(jìn)行了時(shí)空化。 再來(lái)，磁感應(yīng)方程是： 2( ) ( )b u b b b bt (1.2) 這個(gè)方程將磁場(chǎng)進(jìn)行了時(shí)空化。 在 1950 年時(shí) Elsasse 引入 Elsasser 新變量 z u b (1.3) 因?yàn)橛辛诉@個(gè)新的變量， Elsasse 對(duì)動(dòng)量方程和磁感應(yīng)方程進(jìn)行了改 寫(xiě)，并且整合成了一個(gè) 方程： * 2 2( ) ( )A t o tz c z z z P v z v z Ft (1.4) 這個(gè)新方程組“顯式”描述了相反傳播方向的阿爾芬波的非線(xiàn)性相互作用。如圖 1-2 所示 圖 1-2 正向與反向傳播的阿爾文波示意圖 6 1995 年 Goldreich & Shirida 提出基于臨界平衡假設(shè)的理論：波動(dòng)的傳播效應(yīng)和相反波動(dòng)相互作用效應(yīng)相當(dāng)。并且得出如下兩個(gè)方程： ( ) ( )tAz V z w z P (1.5) ( ) ( )tAw V w z w P (1.6) 方程中的 z 和 w 與公式（）中的 z+ 和 z- 一致 VA 均表示是阿爾文速度。 當(dāng)線(xiàn)性傳播項(xiàng)和非線(xiàn)性對(duì)流項(xiàng)的相互作用相當(dāng)?shù)臅r(shí)候，意味著下式： / /AV l b (1.7) 其中 VA 是阿爾文速度， l 是平行磁場(chǎng)的特征尺度；是垂直磁場(chǎng)的特征尺度；并且要求： l (1.8) 式（ 1.7）和（ 1.8）就是所謂的 Critical Balance 也就是臨界平衡。 加上 能量串級(jí)率不隨尺度變化 (為一常數(shù) )，可得： 3V k c o n s t (1.9) AV k V k (1.10) 2( ) /E k V k (1.11) 由以上的條件我們便得到： (1.12) (1.13) 由這兩個(gè)方程我們得到功率譜各向異性的結(jié)論 8,我們基于臨界平衡理論可以做出如下的 PSD(k ,k )圖 : 2()E k k 5 / 3()E k k 7 圖 1-3 基于臨界平衡理論的 PSD(k ,k ) 對(duì)圖 1-3 分別沿著 k積分 ,沿著 k 積分，我們得到圖 1-4 的譜線(xiàn)圖： 圖 1-4 0 度和 90 度的 PSD 圖 8 1.5太陽(yáng)風(fēng)湍流功率譜的各向異性 我們希望通過(guò)觀測(cè)磁場(chǎng)和日心徑向方向不同夾角的功率譜來(lái)觀測(cè)功率譜的各向異性 9，我們先來(lái)看看 Helios 飛船在 0.3AU 的探測(cè)。數(shù)據(jù)經(jīng)過(guò)處理后得到圖 1-5： 圖 1-5 基于 Helios 飛船在 0.3AU 的探測(cè)的 PSD 背景磁場(chǎng)的定義是對(duì)一定時(shí)間取平均值既： B = ；但是這種定義與尺度無(wú)關(guān)，于是我們采用一種新的定義，這種定義使得尺度因子包括進(jìn)背景磁場(chǎng)中，定義式為： (1.14) 我們用 vb 表示觀測(cè)方向與背景磁場(chǎng)的夾角。 PSD(k|) 表示 vb = 0 時(shí)的 PSD 譜線(xiàn)；PSD(k ) 表示 vb = 90 時(shí)的 PSD 譜線(xiàn)； 由圖 1-5 我們可以看出 PSD(k, vb=0)PSD(k, vb=90)和 PSD(k|)PSD(k )，我們對(duì)于這兩組數(shù)據(jù)做出下圖，可以明顯看出這些結(jié)論 圖中紅色線(xiàn)代表 PSD(k )；藍(lán)色線(xiàn)代表 PSD(k|)。 9 圖 1-6 Helios 飛船在 0.3AU 觀測(cè)的 PSD(k|) 和 PSD(k ) 除此之外 ,我們?cè)俳o出 Helios 和 STEREO 飛船在 2007 年 4 月 28 日和 2008 年 2 月 13日兩天觀測(cè)結(jié)果，圖 1-6 表示徑向磁場(chǎng)的觀測(cè)數(shù)據(jù)；中圖表示對(duì)上圖進(jìn)行小波變換后得到的相應(yīng)的功率譜（ PSD）；圖 1-7 是根據(jù)磁場(chǎng)方向與徑向的夾角進(jìn)行按順序重新排列后得到的圖，從圖中很容易看出來(lái)各向異性的特征 10。 10 圖 1-7 Helios、 STEREO 飛船在 2007.04.28 和 2008.02.13 的觀測(cè)結(jié)果圖 這些觀測(cè)都支持了前面所述的臨界平衡假設(shè)的各向異性理論，但是這我們的這篇文章中我們提出這樣的兩個(gè)問(wèn)題；一是這樣的各向異性是一種普適的磁層現(xiàn)象嗎？二是阿爾文波速度場(chǎng)的功率譜是如何的，明顯存在阿爾文波的時(shí)候磁場(chǎng)的功率譜如何。之后我們就來(lái)考察一下這兩個(gè)問(wèn)題。我們采用的是 WIND 飛船在 1995 年 1 月 31 日的數(shù)據(jù)。那天的觀測(cè)證實(shí)太陽(yáng)風(fēng)中明顯存在阿爾文波。 1.6 衛(wèi)星和儀器介紹 WIND11 衛(wèi)星 (圖 1-8；1-9)是長(zhǎng)期觀測(cè)和研究太陽(yáng)風(fēng)的實(shí)驗(yàn)室，它是一個(gè)自旋穩(wěn)定的衛(wèi)星 ,發(fā)射于 1994 年 12 月 .處于日暈軌道繞著第一拉格朗日點(diǎn)運(yùn)動(dòng) ,目的是觀測(cè)穩(wěn)定的太陽(yáng)風(fēng)對(duì)于地球磁層的影響 .wind 和 Geotail, Polar, SoHO and Cluster 一起構(gòu)成了一個(gè)合作學(xué)術(shù)性的用于國(guó)際日地物理項(xiàng)目的衛(wèi)星計(jì)劃 .目的是提高對(duì)于日地間物理關(guān)系的理解。 圖 1-8 WIND 衛(wèi)星介紹 11 Wind 的主要任務(wù)是： 1.提供完整的等離子體 ,能量粒子 ,磁層的磁場(chǎng)和電離層的研究； 2考察出現(xiàn)在近地的太陽(yáng)風(fēng)基本等離子體的過(guò)程； 3 為其他觀測(cè) 提供 1AU 的基線(xiàn)。 圖 1-9 WIND 飛船 我們這篇文章中主要用到了 WIND 上的兩個(gè)儀器，一個(gè)是磁場(chǎng)探測(cè)器 12MFI(Magnetic Field Investigation 圖 1-10；1-11)，另一個(gè)是 3D 等離子分析器 133DP（ 3D Plasma Analyzer 圖 1-12）。 磁場(chǎng)探測(cè)器 MFI 是由兩個(gè)通門(mén)磁力儀構(gòu)成的。其中一個(gè)安裝在飛船上，另一個(gè)安裝在一個(gè)長(zhǎng)達(dá) 12 米的探桿上。這個(gè)儀器測(cè)量直流磁場(chǎng)矢量 ,時(shí)間分辨率為 22 Hz 或者 11Hz，這取決于衛(wèi)星的遙測(cè)模式。 圖 1-10 MFI 的雙單元處理器 圖 1-11 MFI 的雙磁通門(mén)磁力傳感器 12 等離子分析器 3DP 這個(gè)儀器是由 6 個(gè)不同的傳感器組成。他們是 2 個(gè)電子傳感器(EESA)2 個(gè)靜電離子分析器 (PESA)組成。這兩個(gè)傳感器有不同的幾何因子和觀測(cè)范圍。覆蓋能量范圍從 3eV 到 30eV。還有一對(duì)固態(tài)望遠(yuǎn)鏡 (SST)用于測(cè)量能量高達(dá) 400keV 的電子和高達(dá) 6MeV 能量的質(zhì)子。 圖 1-12 Wind 飛船上的 3D 等離子分析器 13 2 數(shù)據(jù)及數(shù)據(jù)處理原理 2.1 數(shù)據(jù)來(lái)源及處理過(guò)程 文章中所用到的數(shù)據(jù)，全部來(lái)自 NASA 官方網(wǎng)站上公布的 1995 年 1 月 31 日的數(shù)據(jù) ,是由 WIND 飛船上的 MFI 和 3DP 收集的。 本文處理數(shù)據(jù)的過(guò)程是（ 1）用 3DP 收集的數(shù)據(jù)繪出速度場(chǎng)的擾動(dòng)，用 MFI 收集的數(shù)據(jù)繪出速度化磁場(chǎng)的擾動(dòng)。將這兩個(gè)擾動(dòng)進(jìn)行比較。（ 2）用 MFI 收集的數(shù)據(jù)繪出局部磁場(chǎng)方向與日心方向的夾角。（ 3）利用磁場(chǎng)擾動(dòng)的數(shù)據(jù)和速度場(chǎng)擾動(dòng)數(shù)據(jù)進(jìn)行小波變換，得到一張 以時(shí)間為橫坐標(biāo)；時(shí)間尺度（ period）為縱坐標(biāo)的功率譜（ 3）將上一步驟得到的兩張功率譜按照時(shí)間相對(duì)應(yīng)的夾角進(jìn)行重新排列，得到一張橫坐標(biāo)為角度的功率譜，縱坐標(biāo)不變。（ 4）通過(guò)重新排列的功率譜就可以進(jìn)行分析并得出結(jié)果。 2.2 小波分析原理 2.2.1 小波分析原理 這部分內(nèi)容參考了 Gtz Paschmann 和 Patrick W. Daly 在 9中的論述。 小波分析在數(shù)據(jù)分析、電磁學(xué)理論、數(shù)據(jù)壓縮等領(lǐng)域都有很重要的應(yīng)用。小波分析的基本觀點(diǎn)是用局限在時(shí)間 t 和頻率 f 的基礎(chǔ)函數(shù)展開(kāi)一個(gè)信號(hào)，因此它們有波包 的性質(zhì)。我們可以選取 Morlet 小波母函數(shù)： 2 0e x p / 2 e x ph t t i t (2.1) 其中 0 是自由參數(shù)。圖 2-1 是由式 2.1 定義的 Morlet 小波母函數(shù)的圖像。 實(shí)際 上， 式 2.1 定義了一族的函數(shù)，稱(chēng)為子函數(shù)。子函數(shù) ()ht 圖 2-1 由式 2.1定義的小波變換母函數(shù) (a)是實(shí)部 (b)是虛部， 0=2 14 1() th t h （ 2.2） 將上式中的 換成頻率 f 1/ ，則式 2.2 變?yōu)椋?22()( ) ( ( ) ) e x p ( ) e x p ( 2 ( ) )2ffth t f h f t f i f t （ 2.3） 0取做 2 。進(jìn)而可以定義對(duì)信號(hào) ()ut 的 Morlet小波變換： *( , ) ( ) ( )fF f u t h t d t （ 2.4） 2.2.2 實(shí)際工作中小波分析的運(yùn)用 正如在 2.2.1中所述，對(duì)于信號(hào) ()ft ，有小波變換： 1 / 2 *( , ) | | ( ) ( )tF s t s f ds （ 2.5） 母函數(shù) ()t 滿(mǎn)足歸一化條件 2| ( ) | 1t d t。且定義一個(gè) C： 2| ( ) |Cd ， ( ) ( ) itt e d t （ 2.6） 并且有： 2 200 /21 / 4 / 2it tt e e e （ 2.7） 在這里常數(shù) 1.06C ，取0 6 。 但在實(shí)際情況下，信號(hào)都是衛(wèi)星測(cè)量的一個(gè)個(gè)數(shù)據(jù)點(diǎn)，不可能做到連續(xù)，因此，實(shí)際小波變換為波變換為： 1 1 / 2 *0()( , ) ( )Nnn k tF s k t s f n t ts （ 2.8） 15 / , 0 , 1 , 2 , , 1t T N n N ， T是記錄的時(shí)間長(zhǎng)度 ,0t CONTOUR,wavelet,time,period ； IDL PLOTS,time,coi,NOCLIP=0 ； ； YPAD = returns the padded time series that was actually used in the 25 ； wavelet transform. ； ； DAUGHTER = if initially set to 1, then return the daughter wavelets. ； This is a complex array of the same size as WAVELET. At each scale ； the daughter wavelet is located in the center of the array. ； ； SIGNIF = output significance levels as a function of PERIOD ； ； FFT_THEOR = output theoretical background spectrum (smoothed by the ； wavelet function), as a function of PERIOD. ； ； ； Defunct INPUTS: ； OCT = the # of octaves to analyze over. ； Largest scale will be S0*2OCT. ； Default is (LOG2(N) - 1). ； VOICE = # of voices in each octave. Default is 8. ； Higher # gives better scale resolution, ； but is slower to plot. ； ； ；- ； ； EXAMPLE: ； ； IDL ntime = 256 ； IDL y = RANDOMN(s,ntime) ；* create a random time series ； IDL dt = 0.25 ； IDL time = FINDGEN(ntime)*dt ；* create the time index ； IDL ； IDL wave = WAVELET(y,dt,PERIOD=period,COI=coi,/PAD,SIGNIF=signif) ； IDL nscale = N_ELEMENTS(period) ； IDL LOADCT,39 ； IDL CONTOUR,ABS(wave)2,time,period, $ ； XSTYLE=1,XTITLE=Time,YTITLE=Period,TITLE=Noise Wavelet, $ ； YRANGE=MAX(period),MIN(period), $ ；* Large-Small period ； /YTYPE, $ ；* make y-axis logarithmic ； NLEVELS=25,/FILL ； IDL ； IDL signif = REBIN(TRANSPOSE(signif),ntime,nscale) ； IDL CONTOUR,ABS(wave)2/signif,time,period, $ ； /OVERPLOT,LEVEL=1.0,C_ANNOT=95% ； IDL PLOTS,time,coi,NOCLIP=0 ；* anything below this line is dubious ； 26 ； ；- ； Copyright (C) 1995-2004, Christopher Torrence and Gilbert P. Compo ； ； This software may be used, copied, or redistributed as long as it is not ； sold and this copyright notice is reproduced on each copy made. ； This routine is provided as is without any express or implied warranties ； whatsoever. ； ； Notice: Please acknowledge the use of the above software in any publications: ； Wavelet software was provided by C. Torrence and G. Compo, ； and is available at URL: /research/wavelets/. ； ； Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to ； Wavelet Analysis. Bull. Amer. Meteor. Soc., 79, 61-78. ； ； Please send a copy of such publications to either C. Torrence or G. Compo: ； Dr. Christopher Torrence Dr. Gilbert P. Compo ； Research Systems, Inc. Climate Diagnostics Center ； 4990 Pearl East Circle 325 Broadway R/CDC1 ； Boulder, CO 80301, USA Boulder, CO 80305-3328, USA ； E-mail: chrisATrsincDOTcom E-mail: compoATcoloradoDOTedu ；- ；- FUNCTION morlet, $ ；* MORLET k0,scale,k,period,coi,dofmin,Cdelta,psi0 IF (k0 EQ -1) THEN k0 = 6d n = N_ELEMENTS(k) expnt = -(scale*k - k0)2/2d*(k GT 0.) dt = 2*!PI/(n*k(1) norm = SQRT(2*!PI*scale/dt)*(!PI(-0.25) ； total energy=N Eqn(7) morlet = norm*EXP(expnt (-100d) morlet = morlet*(expnt GT -100) ； avoid underflow errors morlet = morlet*(k GT 0.) ； Heaviside step function (Morlet is complex) fourier_factor = (4*!PI)/(k0 + SQRT(2+k02) ； Scale-Fourier Sec.3h period = scale*fourier_factor coi = fourier_factor/SQRT(2) ； Cone-of-influence Sec.3g dofmin = 2 ； Degrees of freedom with no smoothing Cdelta = -1 IF (k0 EQ 6) THEN Cdelta = 0.776 ； reconstruction factor psi0 = !PI(-0.25) ； PRINT,scale,n,SQRT(TOTAL(ABS(morlet)2,/DOUBLE) 27 RETURN,morlet END FUNCTION paul, $ ；* PAUL m,scale,k,period,coi,dofmin,Cdelta,psi0 IF (m EQ -1) THEN m = 4d n = N_ELEMENTS(k) expnt = -(scale*k)*(k GT 0.) dt = 2d*!PI/(n*k(1) norm = SQRT(2*!PI*scale/dt)*(2m/SQRT(m*FACTORIAL(2*m-1) paul = norm*(scale*k)m)*EXP(expnt (-100d)*(expnt GT -100) paul = paul*(k GT 0.) fourier_factor = 4*!PI/(2*m+1) period = scale*fourier_factor coi = fourier_factor*SQRT(2) dofmin = 2 ； Degrees of freedom with no smoothing Cdelta = -1 IF (m EQ 4) THEN Cdelta = 1.132 ； reconstruction factor psi0 = 2.m*FACTORIAL(m)/SQRT(!PI*FACTORIAL(2*m) ； PRINT,scale,n,norm,SQRT(TOTAL(paul2,/DOUBLE)*SQRT(n) RETURN,paul END FUNCTION dog, $ ；* DOG m,scale,k,period,coi,dofmin,Cdelta,psi0 IF (m EQ -1) THEN m = 2 n = N_ELEMENTS(k) expnt = -(scale*k)2/2d dt = 2d*!PI/(n*k(1) norm = SQRT(2*!PI*scale/dt)*SQRT(1d/GAMMA(m+0.5) I = DCOMPLEX(0,1) gauss = -norm*(Im)*(scale*k)m*EXP(expnt (-100d)*(expnt GT -100) fourier_factor = 2*!PI*SQRT(2./(2*m+1) period = scale*fourier_factor coi = fourier_factor/SQRT(2) dofmin = 1 ； Degrees of freedom with no smoothing Cdelta = -1 psi0 = -1 IF (m EQ 2) THEN BEGIN Cdelta = 3.541 ； reconstruction factor psi0 = 0.867325 ENDIF 28 IF (m EQ 6) THEN BEGIN Cdelta = 1.966 ； reconstruction factor psi0 = 0.88406 ENDIF ； PRINT,scale,n,norm,SQRT(TOTAL(ABS(gauss)2,/DOUBLE)*SQRT(n) RETURN,gauss END ；* WAVELET FUNCTION wavelet,y1,dt, $ ；* required inputs S0=s0,DJ=dj,J=j, $ ；* optional inputs PAD=pad,MOTHER=mother,PARAM=param, $ VERBOSE=verbose,NO_WAVE=no_wave,RECON=recon, $ LAG1=lag1,SIGLVL=siglvl,DOF=dof,GLOBAL=global, $ ；* optional inputs SCALE=scale,PERIOD=period,YPAD=ypad, $ ；* optional outputs DAUGHTER=daughter,COI=coi, $ SIGNIF=signif,FFT_THEOR=fft_theor, $ OCT=oct,VOICE=voice ；* defunct inputs ON_ERROR,2 r = CHECK_MATH(0,1) n = N_ELEMENTS(y1) n1 = n base2 = FIX(ALOG(n)/ALOG(2) + 0.4999) ； power of 2 nearest to N ；.check keywords & optional inputs IF (N_ELEMENTS(s0) LT 1) THEN s0 = 2.0*dt IF (N_ELEMENTS(voice) EQ 1) THEN dj = 1./voice IF (N_ELEMENTS(dj) LT 1) THEN dj = 1./8 IF (N_ELEMENTS(oct) EQ 1) THEN J = FLOAT(oct)/dj IF (N_ELEMENTS(J) LT 1) THEN J=FIX(ALOG(FLOAT(n)*dt/s0)/ALOG(2)/dj) ；Eqn(10) IF (N_ELEMENTS(mother) LT 1) THEN mother = MORLET IF (N_ELEMENTS(param) LT 1) THEN param = -1 IF (N_ELEMENTS(siglvl) LT 1) THEN siglvl = 0.95 IF (N_ELEMENTS(lag1) LT 1) THEN lag1 = 0.0 lag1 = lag1(0) verbose = KEYWORD_SET(verbose) do_daughter = KEYWORD_SET(daughter) do_wave = NOT KEYWORD_SET(no_wave) recon = KEYWORD_SET(recon) IF KEYWORD_SET(global) THEN MESSAGE, $ Please use WAVE_SIGNIF for global significance tests 29 ；.construct time series to analyze, pad if necessary ypad = y1 - TOTAL(y1)/n ； remove mean IF KEYWORD_SET(pad) THEN BEGIN ； pad with extra zeroes, up to power of 2 ypad = ypad,FLTARR(2L(base2 + 1) - n) n = N_ELEMENTS(ypad) ENDIF ；.construct SCALE array & empty PERIOD & WAVE arrays na = J + 1 ； # of scales scale = DINDGEN(na)*dj ； array of j-values scale = 2d0(scale)*s0 ； array of scales 2j Eqn(9) period = FLTARR(na,/NOZERO) ； empty period array (filled in below) wave = COMPLEXARR(n,na,/NOZERO) ； empty wavelet array IF (do_daughter) THEN daughter = wave ； empty daughter array ；.construct wavenumber array used in transform Eqn(5) k = (DINDGEN(n/2) + 1)*(2*!PI)/(DOUBLE(n)*dt) k = 0d,k,-REVERSE(k(0:(n-1)/2 - 1) ；.compute FFT of the (padded) time series yfft = FFT(ypad,-1,/DOUBLE) ； Eqn(3) IF (verbose) THEN BEGIN ；verbose PRINT PRINT,mother PRINT,#points=,n1, s0=,s0, dj=,dj, J=,FIX(J) IF (n1 NE n) THEN PRINT,(padded with ,n-n1, zeroes) PRINT,j,scale,period,variance,mathflag, $ FORMAT=(/,A3,3A11,A10) ENDIF ；verbose IF (N_ELEMENTS(fft_theor) EQ n) THEN fft_theor_k = fft_theor ELSE $ fft_theor_k = (1-lag12)/(1-2*lag1*COS(k*dt)+lag12) ； Eqn(16) fft_theor = FLTARR(na) ；.loop thru each SCALE FOR a1=0,na-1 DO BEGIN ；scale psi_fft=CALL_FUNCTION(mother, $ param,scale(a1),k,period1,coi,dofmin,Cdelta,psi0) IF (do_wave) THEN $ wave(*,a1) = FFT(yfft*psi_fft,1,/DOUBLE) ；wavelet transformEqn(4) period(a1) = period1 ； save period fft_theor(a1) = TOTAL(ABS(psi_fft)2)*fft_theor_k)/n IF (do_daughter) THEN $ 30 daughter(*,a1) = FFT(psi_fft,1,/DOUBLE) ； save daughter IF (verbose) THEN PRINT,a1,scale(a1),period(a1), $ TOTAL(ABS(wave(*,a1)2),CHECK_MATH(0), $ FORMAT=(I3,3F11.3,I6) ENDFOR ；scale coi = coi*FINDGEN(n1+1)/2),REVERSE(FINDGEN(n1/2)*dt ； COI Sec.3g IF (do_daughter) THEN $ ； shift so DAUGHTERs are in middle of array daughter = daughter(n-n1/2:*,*),daughter(0:n1/2-1,*) ；.significance levels Sec.4 sdev = (MOMENT(y1)(1) fft_theor = sdev*fft_theor ； include time-series variance dof = dofmin signif = fft_theor*CHISQR_CVF(1. - siglvl,dof)/dof ； Eqn(18) IF (recon) THEN BEGIN ； Reconstruction Eqn(11) IF (Cdelta EQ -1) THEN BEGIN y1 = -1 MESSAGE,/INFO, $ Cdelta undefined, cannot reconstruct with this wavelet ENDIF ELSE BEGIN y1=dj*SQRT(dt)/(Cdelta*psi0)*(FLOAT(wave) # (1./SQRT(scale) y1 = y10:n1-1 ENDELSE ENDIF RETURN,wave(0:n1-1,*) ； get rid of padding before returning END 31 附 錄 二 單位轉(zhuǎn)換程序 ；Pro get_BVector_in_AlfvenUnit_vect_WIND_19950131 ； sub_dir_date = 1995-01-31 WIND_Data_Dir = WIND_Data_Dir= WIND_Figure_Dir = WIND_Figure_Dir= SetEnv,WIND_Data_Dir SetEnv,WIND_Figure_Dir Step1: ；= ；Step1: ；- dir_restore = GetEnv(WIND_Data_Dir)+sub_dir_date+ file_restore= B_RTN_3s_arr(time=*).sav file_array = File_Search(dir_restore+file_restore, count=num_files) For i_file=0,num_files-1 Do Begin Print, i_file, file: , i_file, : , file_array(i_file) EndFor i_select = 0 Read, i_select: , i_select file_restore = file_array(i_select) Restore, file_restore, /Verbose ；data_descrip= got from convert_B_GSE_to_B_RTN_WIND_200107.pro ；Save, FileName=dir_save+file_save, $ ； data_descrip, $ ； JulDay_3s_vect, Bxyz_GSE_3s_arr, B_RTN_3s_arr ；- JulDay_vect_B = JulDay_3s_vect ；- dir_restore = GetEnv(WIND_Data_Dir)+sub_dir_date；+MFI file_restore= V_RTN_3s_arr(time=*).sav file_array = File_Search(dir_restore+file_restore, count=num_files) 32 For i_file=0,num_files-1 Do Begin Print, i_file, file: , i_file, : , file_array(i_file) EndFor i_select = 0 Read, i_select: , i_select file_restore = file_array(i_select) Restore, file_restore, /Verbose ；data_descrip= got from convert_V_GSE_to_V_RTN_WIND_200107.pro ；Save, FileName=dir_save+file_save, $ ； data_descrip, $ ； JulDay_3s_vect, Vxyz_GSE_3s_arr, V_RTN_3s_arr ；- JulDay_vect_V = JulDay_3s_vect strpos_tmp = StrPos(file_restore, (time=) strpos_tmp_v2 = StrPos(file_restore, .sav) time_str = StrMid(file_restore, strpos_tmp, strpos_tmp_v2-strpos_tmp) ；- dir_restore = GetEnv(WIND_Data_Dir)+sub_dir_date+3DP file_restore= wi_pm*_3dp_*_*.sav file_array = File_Search(dir_restore+file_restore, count=num_files) For i_file=0,num_files-1 Do Begin Print, i_file, file: , i_file, : , file_array(i_file) EndFor i_select = 0 Read, i_select: , i_select file_restore= file_array(i_select) Restore, file_restore ；data_descrip= got from Read_pm_3dp_WIND_CDF_19950131.pro ；Save, FileName=dir_save+file_save, $ ； data_descrip, $ ； JulDay_vect, Vxyz_GSE_arr, $ ； NumDens_vect, Temper_vect Step2: ；= ；Step2: ；-get new B/V_R/T/N_vect_new ；- JulDay_min = Max(Min(JulDay_vect_B),Min(JulDay_vect_V) JulDay_max = Min(Max(JulDay_vect_B),Max(JulDay_vect_V) 33 num_times_B = N_Elements(JulDay_vect_B) num_times_V = N_Elements(JulDay_vect_V) dJulDay_pix_B = Mean(JulDay_vect_B(1:num_times_B-1)-JulDay_vect_B(0:num_times_B-2) dJulDay_pix_V = Mean(JulDay_vect_V(1:num_times_V-1)-JulDay_vect_V(0:num_times_V-2) dJulDay_pix = Max(dJulDay_pix_B, dJulDay_pix_V) num_times = Floor(JulDay_max-JulDay_min)/dJulDay_pix)+1 ；- JulDay_vect = JulDay_min + Dindgen(num_times)*dJulDay_pix B_RTN_3s_arr= Transpose(B_RTN_3s_arr) V_RTN_3s_arr= Transpose(V_RTN_3s_arr) B_R_vect_old= B_RTN_3s_arr(*,0) B_T_vect_old= B_RTN_3s_arr(*,1) B_N_vect_old= B_RTN_3s_arr(*,2) V_R_vect_old= V_RTN_3s_arr(*,0) V_T_vect_old= V_RTN_3s_arr(*,1) V_N_vect_old= V_RTN_3s_arr(*,2) ；V_R_vect_old= Reform(Vxyz_GSE_3s_arr(0,*) ；V_T_vect_old= Reform(Vxyz_GSE_3s_arr(1,*) ；V_N_vect_old= Reform(Vxyz_GSE_3s_arr(2,*) NumDens_vect_old = NumDens_vect B_R_vect_new= Interpol(B_R_vect_old, JulDay_vect_B, JulDay_vect) B_T_vect_new= Interpol(B_T_vect_old, JulDay_vect_B, JulDay_vect) B_N_vect_new= Interpol(B_N_vect_old, JulDay_vect_B, JulDay_vect) V_R_vect_new= Interpol(V_R_vect_old, JulDay_vect_V, JulDay_vect) V_T_vect_new= Interpol(V_T_vect_old, JulDay_vect_V, JulDay_vect) V_N_vect_new= Interpol(V_N_vect_old, JulDay_vect_V, JulDay_vect) NumDens_vect_new = Interpol(NumDens_vect_old, JulDay_vect_V, JulDay_vect) Step3: ；= ；Step3: ；-get new B_R/T/N_vect_AlfUnit ；- is_instant_or_constant_n = 2 Print, use instant or constant density (1 or 2): , is_instant_or_constant_n is_continue = Read, is_continue: , is_continue ；- 34 If (is_instant_or_constant_n eq 2) Then Begin aver_NumDens = Mean(NumDens_vect_new) NumDens_vect_new_v2 = aver_NumDens + Fltarr(num_times) EndIf If (is_instant_or_constant_n eq 1) Then Begin NumDens_vect_new_v2 = NumDens_vect_new EndIf ；- B_R_vect_AlfUnit = get_AlfvenSpeed(B_R_vect_new, NumDens_vect_new_v2) B_T_vect_AlfUnit = get_AlfvenSpeed(B_T_vect_new, NumDens_vect_new_v2) B_N_vect_AlfUnit = get_AlfvenSpeed(B_N_vect_new, NumDens_vect_new_v2) VA_R_vect_new = B_R_vect_AlfUnit VA_T_vect_new = B_T_vect_AlfUnit VA_N_vect_new = B_N_vect_AlfUnit ；- dir_save = GetEnv(WIND_Data_Dir)+sub_dir_date file_save = V&VA_RTN_3s+time_str+.sav data_descrip= got from get_BVector_in_AlfvenUnit_vect_WIND_19950131.pro Save, FileName=dir_save+file_save, $ data_descrip, $ JulDay_vect, $ B_R_vect_new, B_T_vect_new, B_N_vect_new, $ V_R_vect_new, V_T_vect_new, V_N_vect_new, $ VA_R_vect_new, VA_T_vect_new, VA_N_vect_new, $ NumDens_vect_new End_Program: End 35 附 錄 三 磁場(chǎng)數(shù)據(jù)讀取程序 ；Pro Read_h0_mfi_WIND_CDF_19950131 sub_dir_date = 1995-01-31 Step1: ；= ；Step1: ；- SetEnv, WIND_Data_Dir= dir_read = GetEnv(WIND_Data_Dir)+sub_dir_date+MFI file_read = wi_h0*_mfi_*_*.cdf file_array = File_Search(dir_read+file_read, count=num_files) For i_file=0,num_files-1 Do Begin Print, i_file, file: , i_file, : , file_array(i_file) EndFor i_select = 0 Read, i_select: , i_select file_read = file_array(i_select) file_read = StrMid(file_read, StrLen(dir_read), StrLen(file_read)-StrLen(dir_read) ；- CD, Current=wk_dir CD, dir_read cdf_id = CDF_Open(file_read) ；- CDF_Control, cdf_id, Variable=Epoch3, Get_Var_Info=Info_Epoch3 ；a num_records = Info_Epoch.MaxAllocRec num_records = Info_Epoch3.MaxRec + 1L CDF_VarGet, cdf_id, Epoch3, Epoch_3s_vect, Rec_Count=num_records CDF_VarGet, cdf_id, B3GSE, Bxyz_GSE_3s_arr, Rec_Count=num_records ；- CDF_Control, cdf_id, Variable=Epoch, Get_Var_Info=Info_Epoch num_records = Info_Epoch.MaxRec + 1L CDF_VarGet, cdf_id, Epoch, Epoch_1min_vect, Rec_Count=num_records CDF_VarGet, cdf_id, BGSE, Bxyz_GSE_1min_arr, Rec_Count=num_records CDF_VarGet, cdf_id, PGSE, xyz_GSE_1min_arr, Rec_Count=num_records ；- 36 CDF_Close, cdf_id CD, wk_dir ；- Epoch_3s_vect = Reform(Epoch_3s_vect) epoch_beg = Epoch_3s_vect(0) CDF_Epoch, epoch_beg, year_beg, mon_beg, day_beg, hour_beg, min_beg, sec_beg, milli_beg, $ /Breakdown_Epoch JulDay_beg = JulDay(mon_beg,day_beg,year_beg,hour_beg,min_beg, sec_beg+milli_beg*1.e-3) JulDay_3s_vect = JulDay_beg+(Epoch_3s_vect-Epoch_beg)/(1.e3*24.*60.*60.) ；- Epoch_1min_vect = Reform(Epoch_1min_vect) epoch_beg = Epoch_1min_vect(0) CDF_Epoch, epoch_beg, year_beg, mon_beg, day_beg, hour_beg, min_beg, sec_beg, milli_beg, $ /Breakdown_Epoch JulDay_beg = JulDay(mon_beg,day_beg,year_beg,hour_beg,min_beg, sec_beg+milli_beg*1.e-3) JulDay_1min_vect= JulDay_beg+(Epoch_1min_vect-Epoch_beg)/(1.e3*24.*60.*60.) Step2: ；= ；St