地球物理学毕业论文

题 目 基于临界平衡假设的 太阳风湍流各向异性的观测分析 英文题目 Observation the solar wind turbulence anisotropy based on the critical balance assumption 摘 要 Elsasser 于 1950 年对磁流体方程引入新变量后 Goldreich 和 Shirida 于 1995年提出临界平衡理论假设。并将理论应用于太阳风湍流中得到功率谱各向异性特征和功率谱的谱指标特征。 本文的主要目的是做出速度场的功率谱，进一步检验临界平衡理论假设的正确性与适用性。文中采用了安装在 WIND 卫星上的磁场探测器 MFI 和三维等离子体探测器 3DP 在 1995 年 1 月 31 日所采集到的总磁场和等离子体速度场数据。并且采用基于小波变换的功 率谱分析，做出速度场功率谱密度函数图和磁场功率谱密度函数图。其中横轴为局部磁场和太阳径向向外方向的夹角；纵轴为时间尺度。通过功率谱密度函数图的分析可以得出功率谱随角度变化的分布特征。 通过计算机的计算结果表明：当天的速度场功率谱显式出各向同性的特征，同时磁场的功率谱也同样显示出了各向同性的特征。这说明临界平衡理论假设并不是普遍适用的。 从原始数据看，由于速度扰动和速度化磁场扰动高度相关，笔者认为这可能是产生各向同性的原因之一。遗憾的是本文只针对 1995 年 1 月 31 日的数据进行了考察，无法同其他明显存在阿尔文波 时间段的功率谱进行比较。同时当功率谱为各向同性时，功率谱谱指标特征也有待进一步探究。 关键词 ： 临界平衡；太阳风湍流；功率谱；小波分析；各向异性； 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 卫星介绍 . . 10 2 数据及数据处理原理 . 1 错误 !未定义书签。 2 1 数据来源及处理过程 . 13 2 2 小波分析数据原理 . 13 2 2 1 小波分析原理 . 13 2 2 2 实际工作中小波分析的运用 . 14 2 2 3 功率谱密度 PSD . 15 2 3 功率谱随角度的分布 . 16 3 数据处理及处理结果 . 错误 !未定义书签。 8 3 1 太阳风中的时间序列 . 错误 !未定义书签。 8 3 2 磁场数据处理 . 错误 !未定义书签。 9 3 3 速度场数据处理 . 20 结 论 . 21 致 谢 . 22 参考文献 . 23 附 录 一 . 24 附 录 二 . 32 附 录 三 . 36 1 前 言 1950年， Elsasser对磁流体方程进行了改进，引入了 Elsasser新变量。新变量将磁流体的磁场和动量结合在了一起，不仅将磁流体的动量方程和磁 感应方程结合在一起，同时新方程“显式”的描述了反向传播的阿尔文波的非线性相互作用。 在此基础上 Goldreich 和 Shirida 于 1995年提出了临界平衡理论。他们假设磁流体中波动的传播效应和反向波动相互作用效应相当。并且将理论应用于太阳风中，得到了太阳风湍流各向异性的结论。理论也给出了垂直于磁场方向和平行于磁场方向的功率谱的谱指数分别是 -5/3和 -2 。 Timothy S. Horbury于 2008年的一篇文章中给出了垂直于磁场方向和平行于磁场方向的功率谱的谱指数的观测值，结果与理论符合的很好。 J. J. Podesta于 2009年的一篇文章中进一步细致的研究了局部磁场方向与太阳径向向外方向不同夹角不同时的功率谱谱指标特征，结果同样与理论符合的很好。在 He&Marsch&Tu&Yao&Tian等人共同发表的一篇文章中，通过Helios、 STEREO飞船的观测结果，发现太阳风湍流功率谱的各向异性特征。 很明显，现有的观测结果普遍支持临界平衡理论。不仅观测到了太阳风功率谱的各向异性特征。同时，功率谱谱指数与角度的关系也与理论吻合的很好。但是，临界平衡这个理论现在还处在一个假设性的阶段，并未完全得到肯定。 而且关于速度场的功率谱的研究工作比较少。所以本文的目的在于做出速度场的功率谱，进一步检验临界平衡理论假设的正确性与适用性。这一点对于一个理论的发展是有帮助的。无论最后得到肯定的还是否定的结论，都不仅帮助我们改善我们现有的理论，而且将进一步加深我们对于太阳风湍流特征的认识。从而更好的理解日地之间的空间物理特性。 2 1 背景介绍 1.1 太阳及太阳大气 太阳 1是太阳系的主导恒星，是维持人类生存活动的主要源泉。其质量大约为302 10 kg，半径是 70万千米。日地直线距离是 1.5亿千米，相当于 215个太阳半径。太阳等离子体主要由氢 (约 90%)、氦 (约 10%)以及碳、氮、氧等其他元素 (约 0.1%)组成 章振大 , 1992。太阳的电磁和粒子辐射是影响日地空间环境的主要因素。太阳的电磁辐射覆盖了电磁波谱中从射线到射电波段的极宽的频率范围。同时，太阳连续不断地向外发出数十万度高温的稀薄的磁化等离子体，即以每秒数百公里的速度向外运行的太阳风。 太阳内部从里往外大致可以分为日核、辐射区、对流层。日核 (日心到 0.25Rs)是产能区，太阳在自身引力作用下，物质向日核聚集，导致极高的温度 (107K)，从而不停地进行着由氢核聚变成氦核的热核反应。日核产生的能量在辐射区 (0.25-0.75 Rs)经过多次吸收、散射和再发射逐渐向外传输。辐射区之上是对流层 (0.75 Rs到太阳表面 )，其中的物质处于剧烈的对流状态。 对流层之上便是我们肉眼可见的太阳表面，称为光球。光球及其上的色球、过渡区和日冕合称为太阳大气。图 1.1 中的两条曲线是根据著名的太阳大气 VAL-3C模型 Vernazza et al., 1981计算得到的太阳大气中温度和粒子数密度随高度变化的曲线。通常将太阳大气中温度极小处作为光球和色 球的分界点。占太阳辐射能量中绝大部分的可见光主要来源于光球。对流层的对流运动反映在光球上，主要表现为米粒和超米粒的形式。 光球上方是色球，其温度从 4200 K上升到 20000 K左右。色球辐射主要呈现出亮的或暗的网络结构，网络之间的区域称为网络内区。在高色球中存在着许多针状的结构，被认为是冷而密的色球物质沿磁力线向日冕喷射所形成的。 过渡区是色球之上温度大致在 104 K和 106 K之间的区域。由图 1-1可以看出，在色球和日冕之间的很窄的过渡区内，太阳大气的温度陡增，密度陡降。经过数十年的研究，人们已经认识到 过渡区远非一个静态的分层结构，而是一个磁场和等离子体结构非常不均匀的动态区域。 日冕是太阳大气最外面的一层稀薄的等离子体。 610 K以上的高温导致日冕气体高度电离。日冕辐射主要由如下几部分组成：自由电子对光球辐射的汤姆逊散射形成的 K冕，扩展日冕 (2.5 Rs以上 )中的尘埃散射光球辐射而形成的 F冕，由离子谱线辐射组成的 E冕，以及主要由轫致辐射所产生的 X冕 (软 X射线 )和 R冕 (射电 )。日冕中的结构有辐射很弱的冕洞、活动区里的冕环、极区外呈羽毛状的极羽、以及赤道附近的盔状冕流等。 3 图 1-1 太阳 大气中的温度和粒子数密度随高度的平均变化特征。引自 Peter 2004。 1.2 太阳风及湍流 太阳风是从太阳上层大气射出的超高速等离子体（带电粒子）流。在不是太阳的情况下，这种带电粒子流也常称为“恒星风”。 在太阳的日冕层的高温（几百万开氏度）下，氢、氦等原子已经被电离成带正电的质子、氦原子核和带负电的自由电子等。这些带电粒子运动速度极快，以致不断有带电的粒子挣脱太阳的引力束缚，射向太阳的外围，形成太阳风。 太阳风的速度一般在 200-800km/s。 一般认为在太阳极小期，从太阳的磁场极地附近吹出 的是高速太阳风，从太阳的磁场赤道附近吹出的是低速太阳风。太阳的磁场的活动性是会变化的，周期大约为 11年。 高温导致日冕气体膨胀，连续不断地向外发射等离子体流，到达数个太阳半径的距离后变成超声速的流动，形成太阳风。太阳风主要由电子和质子组成，另有少量粒子和极少量重离子。太阳风把太阳磁场带到行星际空间，形成螺旋状的行星际磁场。 Parker 1958预言太阳风后不久，前苏联和美国的空间飞船探测即证实了太阳风的存在。根据地球轨道 (1 AU)和行星际空间实地测量的流速，可将太阳风分为高速流和低速流两种。 太阳风 一词是在 1950年代被尤金派克提出。但是直到 1960年代才证实了它的存在。长期观测发现，当太阳存在冕洞时，地球附近就能观测到高速的太阳风。因此天文学家认为高速太阳风的产生与冕洞有密切的关系。太阳表面的磁场及等离子体活动对地球有很重要的影响。当太阳发生强烈的活动时，大量的带电粒子随着太阳风吹向地球的两极，就会在两极的电离层引发美丽的极光。 表 1-1列举了 1 AU处观测到的高速流和低速流的平均特性。 4 表 1-1 1 AU处观测到的太阳风高速流和低速流的平均特性。 引自 Axford and McKenzie1997， Xia 2003和 Holzer 2005。 湍流 2，也称为紊流，是流体的一种流动状态。当流速很小时，流体分层流动，互不混合，称为层流，或称为片流；逐渐增加流速，流体的流线开始出现波浪状的摆动，摆动的频率及振幅随流速的增加而增加，此种流况称为过渡流；当流速增加到很大时，流线不再清楚可辨，流场中有许多小漩涡，称为湍流，又称为乱流、扰流或紊流。 而湍流在太阳上和日球层是很普遍的现象。太阳风也是动荡的，并不是平静的流动。可以利用傅里叶谱分析或是小波谱分析来研究太阳风中的湍流 3，来计算其磁场的功率谱 4，结果与早已充分发展的磁流体湍流的结果类似。特别是，功率谱的谱指数接近 -5/35。 1.3 阿尔文波 阿尔文波 6，又称剪切阿尔文波，是等离子体中的一种沿磁场方向传播的波，这种波的频率远低于等离子体的回旋频率，是一种线偏振的低频横波。处在磁场中的导电流体在垂直于磁场的方向上受到局部扰动时，沿着磁感线方向的磁张力提供恢复力，就会激发阿尔文波。阿尔文波是由瑞典物理学家汉尼斯阿尔文首先预言的，因此得名。后来隆德奎斯特使用 1特斯拉左右的磁场在水银中观察到了阿尔文波，列纳尔特使用 液态钠也证实了阿尔文波的存在。阿尔文波的色散关系为： 这样的波称为斜阿尔文波。 =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 均表示是阿尔文速度。 当线性传播项和非线性对流项的相互作用相当的时候，意味着下式： / /AV l b (1.7) 其中 VA 是阿尔文速度， l 是平行磁场的特征尺度；是垂直磁场的特征尺度；并且要求： l (1.8) 式（ 1.7）和（ 1.8）就是所谓的 Critical Balance 也就是临界平衡。 加上 能量串级率不随尺度变化 (为一常数 )，可得： 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 的探测。数据经过处理后得到图 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 )，我们对于这两组数据做出下图，可以明显看出这些结论 图中红色线代表 PSD(k )；蓝色线代表 PSD(k|)。 9 图 1-6 Helios 飞船在 0.3AU 观测的 PSD(k|) 和 PSD(k ) 除此之外 ,我们再给出 Helios 和 STEREO 飞船在 2007 年 4 月 28 日和 2008 年 2 月 13日两天观测结果，图 1-6 表示径向磁场的观测数据；中图表示对上图进行小波变换后得到的相应的功率谱（ PSD）；图 1-7 是根据磁场方向与径向的夹角进行按顺序重新排列后得到的图，从图中很容易看出来各向异性的特征 10。 10 图 1-7 Helios、 STEREO 飞船在 2007.04.28 和 2008.02.13 的观测结果图 这些观测都支持了前面所述的临界平衡假设的各向异性理论，但是这我们的这篇文章中我们提出这样的两个问题；一是这样的各向异性是一种普适的磁层现象吗？二是阿尔文波速度场的功率谱是如何的，明显存在阿尔文波的时候磁场的功率谱如何。之后我们就来考察一下这两个问题。我们采用的是 WIND 飞船在 1995 年 1 月 31 日的数据。那天的观测证实太阳风中明显存在阿尔文波。 1.6 卫星和仪器介绍 WIND11 卫星 (图 1-8；1-9)是长期观测和研究太阳风的实验室，它是一个自旋稳定的卫星 ,发射于 1994 年 12 月 .处于日晕轨道绕着第一拉格朗日点运动 ,目的是观测稳定的太阳风对于地球磁层的影响 .wind 和 Geotail, Polar, SoHO and Cluster 一起构成了一个合作学术性的用于国际日地物理项目的卫星计划 .目的是提高对于日地间物理关系的理解。 图 1-8 WIND 卫星介绍 11 Wind 的主要任务是： 1.提供完整的等离子体 ,能量粒子 ,磁层的磁场和电离层的研究； 2考察出现在近地的太阳风基本等离子体的过程； 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，这取决于卫星的遥测模式。 图 1-10 MFI 的双单元处理器 图 1-11 MFI 的双磁通门磁力传感器 12 等离子分析器 3DP 这个仪器是由 6 个不同的传感器组成。他们是 2 个电子传感器(EESA)2 个静电离子分析器 (PESA)组成。这两个传感器有不同的几何因子和观测范围。覆盖能量范围从 3eV 到 30eV。还有一对固态望远镜 (SST)用于测量能量高达 400keV 的电子和高达 6MeV 能量的质子。 图 1-12 Wind 飞船上的 3D 等离子分析器 13 2 数据及数据处理原理 2.1 数据来源及处理过程 文章中所用到的数据，全部来自 NASA 官方网站上公布的 1995 年 1 月 31 日的数据 ,是由 WIND 飞船上的 MFI 和 3DP 收集的。 本文处理数据的过程是（ 1）用 3DP 收集的数据绘出速度场的扰动，用 MFI 收集的数据绘出速度化磁场的扰动。将这两个扰动进行比较。（ 2）用 MFI 收集的数据绘出局部磁场方向与日心方向的夹角。（ 3）利用磁场扰动的数据和速度场扰动数据进行小波变换，得到一张 以时间为横坐标；时间尺度（ period）为纵坐标的功率谱（ 3）将上一步骤得到的两张功率谱按照时间相对应的夹角进行重新排列，得到一张横坐标为角度的功率谱，纵坐标不变。（ 4）通过重新排列的功率谱就可以进行分析并得出结果。 2.2 小波分析原理 2.2.1 小波分析原理 这部分内容参考了 Gtz Paschmann 和 Patrick W. Daly 在 9中的论述。 小波分析在数据分析、电磁学理论、数据压缩等领域都有很重要的应用。小波分析的基本观点是用局限在时间 t 和频率 f 的基础函数展开一个信号，因此它们有波包 的性质。我们可以选取 Morlet 小波母函数： 2 0e x p / 2 e x ph t t i t (2.1) 其中 0 是自由参数。图 2-1 是由式 2.1 定义的 Morlet 小波母函数的图像。 实际 上， 式 2.1 定义了一族的函数，称为子函数。子函数 ()ht 图 2-1 由式 2.1定义的小波变换母函数 (a)是实部 (b)是虚部， 0=2 14 1() th t h （ 2.2） 将上式中的 换成频率 f 1/ ，则式 2.2 变为： 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） 母函数 ()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） 在这里常数 1.06C ，取0 6 。 但在实际情况下，信号都是卫星测量的一个个数据点，不可能做到连续，因此，实际小波变换为波变换为： 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 附 录 三 磁场数据读取程序 ；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