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

題 目 基于臨界平衡假設的 太陽風湍流各向異性的觀測分析 英文題目 Observation the solar wind turbulence anisotropy based on the critical balance assumption 摘 要 Elsasser 于 1950 年對磁流體方程引入新變量后 Goldreich 和 Shirida 于 1995年提出臨界平衡理論假設。并將理論應用于太陽風湍流中得到功率譜各向異性特征和功率譜的譜指標特征。 本文的主要目的是做出速度場的功率譜，進一步檢驗臨界平衡理論假設的正確性與適用性。文中采用了安裝在 WIND 衛(wèi)星上的磁場探測器 MFI 和三維等離子體探測器 3DP 在 1995 年 1 月 31 日所采集到的總磁場和等離子體速度場數(shù)據(jù)。并且采用基于小波變換的功 率譜分析，做出速度場功率譜密度函數(shù)圖和磁場功率譜密度函數(shù)圖。其中橫軸為局部磁場和太陽徑向向外方向的夾角；縱軸為時間尺度。通過功率譜密度函數(shù)圖的分析可以得出功率譜隨角度變化的分布特征。 通過計算機的計算結果表明：當天的速度場功率譜顯式出各向同性的特征，同時磁場的功率譜也同樣顯示出了各向同性的特征。這說明臨界平衡理論假設并不是普遍適用的。 從原始數(shù)據(jù)看，由于速度擾動和速度化磁場擾動高度相關，筆者認為這可能是產(chǎn)生各向同性的原因之一。遺憾的是本文只針對 1995 年 1 月 31 日的數(shù)據(jù)進行了考察，無法同其他明顯存在阿爾文波 時間段的功率譜進行比較。同時當功率譜為各向同性時，功率譜譜指標特征也有待進一步探究。 關鍵詞 ： 臨界平衡；太陽風湍流；功率譜；小波分析；各向異性； 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； 目 錄 前 言 . 錯誤 !未定義書簽。 1 背景介紹 . 錯誤 !未定義書簽。 1 1 太陽大氣與太陽風 . 錯誤 !未定義書簽。 1 2 太陽風及湍流 . 3 1 3 阿爾文波 . 4 1 4 臨界平衡 . 5 1 5 太陽風湍流功率譜的各向異性 . 8 1 6 衛(wèi)星介紹 . . 10 2 數(shù)據(jù)及數(shù)據(jù)處理原理 . 1 錯誤 !未定義書簽。 2 1 數(shù)據(jù)來源及處理過程 . 13 2 2 小波分析數(shù)據(jù)原理 . 13 2 2 1 小波分析原理 . 13 2 2 2 實際工作中小波分析的運用 . 14 2 2 3 功率譜密度 PSD . 15 2 3 功率譜隨角度的分布 . 16 3 數(shù)據(jù)處理及處理結果 . 錯誤 !未定義書簽。 8 3 1 太陽風中的時間序列 . 錯誤 !未定義書簽。 8 3 2 磁場數(shù)據(jù)處理 . 錯誤 !未定義書簽。 9 3 3 速度場數(shù)據(jù)處理 . 20 結 論 . 21 致 謝 . 22 參考文獻 . 23 附 錄 一 . 24 附 錄 二 . 32 附 錄 三 . 36 1 前 言 1950年， Elsasser對磁流體方程進行了改進，引入了 Elsasser新變量。新變量將磁流體的磁場和動量結合在了一起，不僅將磁流體的動量方程和磁 感應方程結合在一起，同時新方程“顯式”的描述了反向傳播的阿爾文波的非線性相互作用。 在此基礎上 Goldreich 和 Shirida 于 1995年提出了臨界平衡理論。他們假設磁流體中波動的傳播效應和反向波動相互作用效應相當。并且將理論應用于太陽風中，得到了太陽風湍流各向異性的結論。理論也給出了垂直于磁場方向和平行于磁場方向的功率譜的譜指數(shù)分別是 -5/3和 -2 。 Timothy S. Horbury于 2008年的一篇文章中給出了垂直于磁場方向和平行于磁場方向的功率譜的譜指數(shù)的觀測值，結果與理論符合的很好。 J. J. Podesta于 2009年的一篇文章中進一步細致的研究了局部磁場方向與太陽徑向向外方向不同夾角不同時的功率譜譜指標特征，結果同樣與理論符合的很好。在 He&Marsch&Tu&Yao&Tian等人共同發(fā)表的一篇文章中，通過Helios、 STEREO飛船的觀測結果，發(fā)現(xiàn)太陽風湍流功率譜的各向異性特征。 很明顯，現(xiàn)有的觀測結果普遍支持臨界平衡理論。不僅觀測到了太陽風功率譜的各向異性特征。同時，功率譜譜指數(shù)與角度的關系也與理論吻合的很好。但是，臨界平衡這個理論現(xiàn)在還處在一個假設性的階段，并未完全得到肯定。 而且關于速度場的功率譜的研究工作比較少。所以本文的目的在于做出速度場的功率譜，進一步檢驗臨界平衡理論假設的正確性與適用性。這一點對于一個理論的發(fā)展是有幫助的。無論最后得到肯定的還是否定的結論，都不僅幫助我們改善我們現(xiàn)有的理論，而且將進一步加深我們對于太陽風湍流特征的認識。從而更好的理解日地之間的空間物理特性。 2 1 背景介紹 1.1 太陽及太陽大氣 太陽 1是太陽系的主導恒星，是維持人類生存活動的主要源泉。其質量大約為302 10 kg，半徑是 70萬千米。日地直線距離是 1.5億千米，相當于 215個太陽半徑。太陽等離子體主要由氫 (約 90%)、氦 (約 10%)以及碳、氮、氧等其他元素 (約 0.1%)組成 章振大 , 1992。太陽的電磁和粒子輻射是影響日地空間環(huán)境的主要因素。太陽的電磁輻射覆蓋了電磁波譜中從射線到射電波段的極寬的頻率范圍。同時，太陽連續(xù)不斷地向外發(fā)出數(shù)十萬度高溫的稀薄的磁化等離子體，即以每秒數(shù)百公里的速度向外運行的太陽風。 太陽內部從里往外大致可以分為日核、輻射區(qū)、對流層。日核 (日心到 0.25Rs)是產(chǎn)能區(qū)，太陽在自身引力作用下，物質向日核聚集，導致極高的溫度 (107K)，從而不停地進行著由氫核聚變成氦核的熱核反應。日核產(chǎn)生的能量在輻射區(qū) (0.25-0.75 Rs)經(jīng)過多次吸收、散射和再發(fā)射逐漸向外傳輸。輻射區(qū)之上是對流層 (0.75 Rs到太陽表面 )，其中的物質處于劇烈的對流狀態(tài)。 對流層之上便是我們肉眼可見的太陽表面，稱為光球。光球及其上的色球、過渡區(qū)和日冕合稱為太陽大氣。圖 1.1 中的兩條曲線是根據(jù)著名的太陽大氣 VAL-3C模型 Vernazza et al., 1981計算得到的太陽大氣中溫度和粒子數(shù)密度隨高度變化的曲線。通常將太陽大氣中溫度極小處作為光球和色 球的分界點。占太陽輻射能量中絕大部分的可見光主要來源于光球。對流層的對流運動反映在光球上，主要表現(xiàn)為米粒和超米粒的形式。 光球上方是色球，其溫度從 4200 K上升到 20000 K左右。色球輻射主要呈現(xiàn)出亮的或暗的網(wǎng)絡結構，網(wǎng)絡之間的區(qū)域稱為網(wǎng)絡內區(qū)。在高色球中存在著許多針狀的結構，被認為是冷而密的色球物質沿磁力線向日冕噴射所形成的。 過渡區(qū)是色球之上溫度大致在 104 K和 106 K之間的區(qū)域。由圖 1-1可以看出，在色球和日冕之間的很窄的過渡區(qū)內，太陽大氣的溫度陡增，密度陡降。經(jīng)過數(shù)十年的研究，人們已經(jīng)認識到 過渡區(qū)遠非一個靜態(tài)的分層結構，而是一個磁場和等離子體結構非常不均勻的動態(tài)區(qū)域。 日冕是太陽大氣最外面的一層稀薄的等離子體。 610 K以上的高溫導致日冕氣體高度電離。日冕輻射主要由如下幾部分組成：自由電子對光球輻射的湯姆遜散射形成的 K冕，擴展日冕 (2.5 Rs以上 )中的塵埃散射光球輻射而形成的 F冕，由離子譜線輻射組成的 E冕，以及主要由軔致輻射所產(chǎn)生的 X冕 (軟 X射線 )和 R冕 (射電 )。日冕中的結構有輻射很弱的冕洞、活動區(qū)里的冕環(huán)、極區(qū)外呈羽毛狀的極羽、以及赤道附近的盔狀冕流等。 3 圖 1-1 太陽 大氣中的溫度和粒子數(shù)密度隨高度的平均變化特征。引自 Peter 2004。 1.2 太陽風及湍流 太陽風是從太陽上層大氣射出的超高速等離子體（帶電粒子）流。在不是太陽的情況下，這種帶電粒子流也常稱為“恒星風”。 在太陽的日冕層的高溫（幾百萬開氏度）下，氫、氦等原子已經(jīng)被電離成帶正電的質子、氦原子核和帶負電的自由電子等。這些帶電粒子運動速度極快，以致不斷有帶電的粒子掙脫太陽的引力束縛，射向太陽的外圍，形成太陽風。 太陽風的速度一般在 200-800km/s。 一般認為在太陽極小期，從太陽的磁場極地附近吹出 的是高速太陽風，從太陽的磁場赤道附近吹出的是低速太陽風。太陽的磁場的活動性是會變化的，周期大約為 11年。 高溫導致日冕氣體膨脹，連續(xù)不斷地向外發(fā)射等離子體流，到達數(shù)個太陽半徑的距離后變成超聲速的流動，形成太陽風。太陽風主要由電子和質子組成，另有少量粒子和極少量重離子。太陽風把太陽磁場帶到行星際空間，形成螺旋狀的行星際磁場。 Parker 1958預言太陽風后不久，前蘇聯(lián)和美國的空間飛船探測即證實了太陽風的存在。根據(jù)地球軌道 (1 AU)和行星際空間實地測量的流速，可將太陽風分為高速流和低速流兩種。 太陽風 一詞是在 1950年代被尤金派克提出。但是直到 1960年代才證實了它的存在。長期觀測發(fā)現(xiàn)，當太陽存在冕洞時，地球附近就能觀測到高速的太陽風。因此天文學家認為高速太陽風的產(chǎn)生與冕洞有密切的關系。太陽表面的磁場及等離子體活動對地球有很重要的影響。當太陽發(fā)生強烈的活動時，大量的帶電粒子隨著太陽風吹向地球的兩極，就會在兩極的電離層引發(fā)美麗的極光。 表 1-1列舉了 1 AU處觀測到的高速流和低速流的平均特性。 4 表 1-1 1 AU處觀測到的太陽風高速流和低速流的平均特性。 引自 Axford and McKenzie1997， Xia 2003和 Holzer 2005。 湍流 2，也稱為紊流，是流體的一種流動狀態(tài)。當流速很小時，流體分層流動，互不混合，稱為層流，或稱為片流；逐漸增加流速，流體的流線開始出現(xiàn)波浪狀的擺動，擺動的頻率及振幅隨流速的增加而增加，此種流況稱為過渡流；當流速增加到很大時，流線不再清楚可辨，流場中有許多小漩渦，稱為湍流，又稱為亂流、擾流或紊流。 而湍流在太陽上和日球層是很普遍的現(xiàn)象。太陽風也是動蕩的，并不是平靜的流動?？梢岳酶道锶~譜分析或是小波譜分析來研究太陽風中的湍流 3，來計算其磁場的功率譜 4，結果與早已充分發(fā)展的磁流體湍流的結果類似。特別是，功率譜的譜指數(shù)接近 -5/35。 1.3 阿爾文波 阿爾文波 6，又稱剪切阿爾文波，是等離子體中的一種沿磁場方向傳播的波，這種波的頻率遠低于等離子體的回旋頻率，是一種線偏振的低頻橫波。處在磁場中的導電流體在垂直于磁場的方向上受到局部擾動時，沿著磁感線方向的磁張力提供恢復力，就會激發(fā)阿爾文波。阿爾文波是由瑞典物理學家漢尼斯阿爾文首先預言的，因此得名。后來隆德奎斯特使用 1特斯拉左右的磁場在水銀中觀察到了阿爾文波，列納爾特使用 液態(tài)鈉也證實了阿爾文波的存在。阿爾文波的色散關系為： 這樣的波稱為斜阿爾文波。 =0 時是沿著磁感線的方向傳播的，稱為阿爾文波。此時阿爾文波的相速度為： 稱為斜傳阿爾文速度，其中 0 是等離子體的磁導率， 0是離子的密度。在垂直于磁感線的方向上阿爾文波不能傳播。 5 1.4 臨界平衡 關于太陽風的傳播， 1995 年 Goldreich 和 Shirida 提出基于臨界平衡假設的理論，這里做一個簡單的介紹。 我們先來看一下非線性簡化磁流體 (MHD)方程的新形式 7。 首先，動量方程是： 2( ) ( )t o tu u u P v u b bt (1.1) 這個方程將速度進行了時空化。 再來，磁感應方程是： 2( ) ( )b u b b b bt (1.2) 這個方程將磁場進行了時空化。 在 1950 年時 Elsasse 引入 Elsasser 新變量 z u b (1.3) 因為有了這個新的變量， Elsasse 對動量方程和磁感應方程進行了改 寫，并且整合成了一個 方程： * 2 2( ) ( )A t o tz c z z z P v z v z Ft (1.4) 這個新方程組“顯式”描述了相反傳播方向的阿爾芬波的非線性相互作用。如圖 1-2 所示 圖 1-2 正向與反向傳播的阿爾文波示意圖 6 1995 年 Goldreich & Shirida 提出基于臨界平衡假設的理論：波動的傳播效應和相反波動相互作用效應相當。并且得出如下兩個方程： ( ) ( )tAz V z w z P (1.5) ( ) ( )tAw V w z w P (1.6) 方程中的 z 和 w 與公式（）中的 z+ 和 z- 一致 VA 均表示是阿爾文速度。 當線性傳播項和非線性對流項的相互作用相當?shù)臅r候，意味著下式： / /AV l b (1.7) 其中 VA 是阿爾文速度， l 是平行磁場的特征尺度；是垂直磁場的特征尺度；并且要求： l (1.8) 式（ 1.7）和（ 1.8）就是所謂的 Critical Balance 也就是臨界平衡。 加上 能量串級率不隨尺度變化 (為一常數(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) 由這兩個方程我們得到功率譜各向異性的結論 8,我們基于臨界平衡理論可以做出如下的 PSD(k ,k )圖 : 2()E k k 5 / 3()E k k 7 圖 1-3 基于臨界平衡理論的 PSD(k ,k ) 對圖 1-3 分別沿著 k積分 ,沿著 k 積分，我們得到圖 1-4 的譜線圖： 圖 1-4 0 度和 90 度的 PSD 圖 8 1.5太陽風湍流功率譜的各向異性 我們希望通過觀測磁場和日心徑向方向不同夾角的功率譜來觀測功率譜的各向異性 9，我們先來看看 Helios 飛船在 0.3AU 的探測。數(shù)據(jù)經(jīng)過處理后得到圖 1-5： 圖 1-5 基于 Helios 飛船在 0.3AU 的探測的 PSD 背景磁場的定義是對一定時間取平均值既： B = ；但是這種定義與尺度無關，于是我們采用一種新的定義，這種定義使得尺度因子包括進背景磁場中，定義式為： (1.14) 我們用 vb 表示觀測方向與背景磁場的夾角。 PSD(k|) 表示 vb = 0 時的 PSD 譜線；PSD(k ) 表示 vb = 90 時的 PSD 譜線； 由圖 1-5 我們可以看出 PSD(k, vb=0)PSD(k, vb=90)和 PSD(k|)PSD(k )，我們對于這兩組數(shù)據(jù)做出下圖，可以明顯看出這些結論 圖中紅色線代表 PSD(k )；藍色線代表 PSD(k|)。 9 圖 1-6 Helios 飛船在 0.3AU 觀測的 PSD(k|) 和 PSD(k ) 除此之外 ,我們再給出 Helios 和 STEREO 飛船在 2007 年 4 月 28 日和 2008 年 2 月 13日兩天觀測結果，圖 1-6 表示徑向磁場的觀測數(shù)據(jù)；中圖表示對上圖進行小波變換后得到的相應的功率譜（ PSD）；圖 1-7 是根據(jù)磁場方向與徑向的夾角進行按順序重新排列后得到的圖，從圖中很容易看出來各向異性的特征 10。 10 圖 1-7 Helios、 STEREO 飛船在 2007.04.28 和 2008.02.13 的觀測結果圖 這些觀測都支持了前面所述的臨界平衡假設的各向異性理論，但是這我們的這篇文章中我們提出這樣的兩個問題；一是這樣的各向異性是一種普適的磁層現(xiàn)象嗎？二是阿爾文波速度場的功率譜是如何的，明顯存在阿爾文波的時候磁場的功率譜如何。之后我們就來考察一下這兩個問題。我們采用的是 WIND 飛船在 1995 年 1 月 31 日的數(shù)據(jù)。那天的觀測證實太陽風中明顯存在阿爾文波。 1.6 衛(wèi)星和儀器介紹 WIND11 衛(wèi)星 (圖 1-8；1-9)是長期觀測和研究太陽風的實驗室，它是一個自旋穩(wěn)定的衛(wèi)星 ,發(fā)射于 1994 年 12 月 .處于日暈軌道繞著第一拉格朗日點運動 ,目的是觀測穩(wěn)定的太陽風對于地球磁層的影響 .wind 和 Geotail, Polar, SoHO and Cluster 一起構成了一個合作學術性的用于國際日地物理項目的衛(wèi)星計劃 .目的是提高對于日地間物理關系的理解。 圖 1-8 WIND 衛(wèi)星介紹 11 Wind 的主要任務是： 1.提供完整的等離子體 ,能量粒子 ,磁層的磁場和電離層的研究； 2考察出現(xiàn)在近地的太陽風基本等離子體的過程； 3 為其他觀測 提供 1AU 的基線。 圖 1-9 WIND 飛船 我們這篇文章中主要用到了 WIND 上的兩個儀器，一個是磁場探測器 12MFI(Magnetic Field Investigation 圖 1-10；1-11)，另一個是 3D 等離子分析器 133DP（ 3D Plasma Analyzer 圖 1-12）。 磁場探測器 MFI 是由兩個通門磁力儀構成的。其中一個安裝在飛船上，另一個安裝在一個長達 12 米的探桿上。這個儀器測量直流磁場矢量 ,時間分辨率為 22 Hz 或者 11Hz，這取決于衛(wèi)星的遙測模式。 圖 1-10 MFI 的雙單元處理器 圖 1-11 MFI 的雙磁通門磁力傳感器 12 等離子分析器 3DP 這個儀器是由 6 個不同的傳感器組成。他們是 2 個電子傳感器(EESA)2 個靜電離子分析器 (PESA)組成。這兩個傳感器有不同的幾何因子和觀測范圍。覆蓋能量范圍從 3eV 到 30eV。還有一對固態(tài)望遠鏡 (SST)用于測量能量高達 400keV 的電子和高達 6MeV 能量的質子。 圖 1-12 Wind 飛船上的 3D 等離子分析器 13 2 數(shù)據(jù)及數(shù)據(jù)處理原理 2.1 數(shù)據(jù)來源及處理過程 文章中所用到的數(shù)據(jù)，全部來自 NASA 官方網(wǎng)站上公布的 1995 年 1 月 31 日的數(shù)據(jù) ,是由 WIND 飛船上的 MFI 和 3DP 收集的。 本文處理數(shù)據(jù)的過程是（ 1）用 3DP 收集的數(shù)據(jù)繪出速度場的擾動，用 MFI 收集的數(shù)據(jù)繪出速度化磁場的擾動。將這兩個擾動進行比較。（ 2）用 MFI 收集的數(shù)據(jù)繪出局部磁場方向與日心方向的夾角。（ 3）利用磁場擾動的數(shù)據(jù)和速度場擾動數(shù)據(jù)進行小波變換，得到一張 以時間為橫坐標；時間尺度（ period）為縱坐標的功率譜（ 3）將上一步驟得到的兩張功率譜按照時間相對應的夾角進行重新排列，得到一張橫坐標為角度的功率譜，縱坐標不變。（ 4）通過重新排列的功率譜就可以進行分析并得出結果。 2.2 小波分析原理 2.2.1 小波分析原理 這部分內容參考了 Gtz Paschmann 和 Patrick W. Daly 在 9中的論述。 小波分析在數(shù)據(jù)分析、電磁學理論、數(shù)據(jù)壓縮等領域都有很重要的應用。小波分析的基本觀點是用局限在時間 t 和頻率 f 的基礎函數(shù)展開一個信號，因此它們有波包 的性質。我們可以選取 Morlet 小波母函數(shù)： 2 0e x p / 2 e x ph t t i t (2.1) 其中 0 是自由參數(shù)。圖 2-1 是由式 2.1 定義的 Morlet 小波母函數(shù)的圖像。 實際 上， 式 2.1 定義了一族的函數(shù)，稱為子函數(shù)。子函數(shù) ()ht 圖 2-1 由式 2.1定義的小波變換母函數(shù) (a)是實部 (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 。進而可以定義對信號 ()ut 的 Morlet小波變換： *( , ) ( ) ( )fF f u t h t d t （ 2.4） 2.2.2 實際工作中小波分析的運用 正如在 2.2.1中所述，對于信號 ()ft ，有小波變換： 1 / 2 *( , ) | | ( ) ( )tF s t s f ds （ 2.5） 母函數(shù) ()t 滿足歸一化條件 2| ( ) | 1t d t。且定義一個 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 。 但在實際情況下，信號都是衛(wèi)星測量的一個個數(shù)據(jù)點，不可能做到連續(xù)，因此，實際小波變換為波變換為： 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是記錄的時間長度 ,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 附 錄 二 單位轉換程序 ；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 附 錄 三 磁場數(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