0 引言

2018年2月2日中国自主研制的首颗地震电磁监测试验卫星——张衡一号(China Seismo-Electromagnetic Satellite，下文简称CSES)在中国酒泉卫星发射中心成功发射入轨.该卫星科学目标是获取全球电磁场、电离层等离子体、高能粒子观测数据，对中国及周边区域开展电离层动态实时监测和地震前兆跟踪监测，探索地震电离层扰动机制(Shen et al., 2018a, b).

CSES卫星运行轨道为507 km高度的圆极地轨道，轨道倾角为97°，重访周期为5天，这个轨道空间属于地球低轨道(Low Earth Orbit)卫星空间.在此低轨道卫星空间内，CSES卫星在高纬度地区可以穿越等离子体层的层顶(Zhima et al., 2013)、极尖区，在中低纬度地区经过粒子沉降区等，而这些区域是岩石圈-大气层-电离层和电离层-磁层圈层耦合的重要作用区域(曹晋滨等，2009).在这个区域内，来自太阳、地磁场的扰动与地震、闪电、人类活动等造成的扰动混合在一起，构成了复杂的扰动体系，产生了各种类型的电磁波动现象(Parrot et al., 2006a; Santolík et al., 2006; Zhima et al., 2014, 2015).

根据DEMETER卫星在顶部电离层长达六年的翔实观测可知顶部电离层区域的ULF/ELF/VLF电磁波动主要分为上行和下行两种类型(Zhima et al., 2013, 2017).上行电磁波动是从岩石圈、大气层向上传播的电磁波动，根据波矢量分析，这些电磁波动的波矢量和Poynting能流从地球方向向外空间传播(L shell增大的方向).下行电磁波动是从等离子体层或内磁层向下传播到顶部电离层的电磁波动(Chen et al., 2017; Santolík et al., 2006; Zhima et al., 2017)，这种类型电磁波动的波矢量和Poynting能流从外空间向地球方向传播(L shell减小的方向).

前人研究结果也证实了这个区域内卫星观测的地震电离层异常(Parrot et al., 2006a; Rozhnoi et al., 2010; Akhoondzadeh et al., 2010; Zhima et al., 2012; Tao et al., 2018；张学民等，2009；于海雁和乔晓林，2009；曾中超等，2009；朱涛和王兰炜，2011；泽仁志玛等，2012).岩石圈发生的强烈地震可能在ULF/ELF/VLF频段上激发不同类型的电磁扰动，如汶川地震(曾中超等，2009)、智利地震(Zhang et al., 2011)、玉树地震(Zhima et al., 2012)前后DEMETER卫星在ELF/VLF频段上观测到了强烈的电磁波动.统计研究也发现强震孕育区上空ELF/VLF频段的变化磁场，在地震发生前后相对于其背景场存在一定的扰动幅度，而在全球范围内随机选择的一组无地震发生地点上空并没有发现类似扰动(泽仁志玛等，2012).然而无强震发生期间，在这个空间内ULF/ELF/VLF频段(3 kHz)以下的电磁辐射非常普遍，即使在地磁场宁静情况下电磁辐射依然很活跃，其空间展布从中低纬度一直到高纬度地区(Zhima et al., 2014).

研究空间的电磁波动传播特征，需要利用电、磁场多分量的波形数据进行频谱分析和波矢量分析.CSES卫星的电磁场观测提供了ULF/ELF/VLF频段六分量电磁波形数据，这为我们计算电磁波的传播参数提供了基本的数据条件.通常对时域的电磁数据我们需要采用频域分析，因为在频域内可以更真实地反映波的信息，其次利用电磁场多分量数据开展波矢量分析，根据电磁理论知识可知电磁波的极化平面和波矢方向相互垂直，通过波矢量分析可获得波的振幅、相位、传播角、频移、极化等参数.但是波矢量分析只能确定波矢量的平行两个方向，存在±180°的方位差，因此通常采用Poynting矢量的能流方向辅助判定，弥补波矢量分析无法确定波的具体传播方向的缺点.本研究主要利用张衡一号卫星记录的电磁场波形数据，开发了相应频谱、波矢分析算法工具，并对CSES卫星记录的几种典型波动事件开展了分析研究.

1 卫星数据简介

本研究主要采用了CSES和DEMETER卫星观测数据，利用DEMETER卫星的主要目的是验证关于波矢分析算法的正确性，以及检验CSES卫星对相同事件的观测能力.CSES卫星的成功发射为开展全球空间电磁场、电离层等离子体、高能粒子沉降等物理现象的监测以及地震机理研究、空间环境监测和地球系统科学研究提供了新的技术手段(Shen et al., 2018a, b).CSES卫星在运行期间以巡查模式为主，经过中国区域或者全球地震带上空则调整为详查模式，日侧和夜侧的升降交点地方时分别为02:00、14:00 LT，每个轨道周期大概为94.6 min，重访周期为5天.CSES搭载了感应式磁力仪、高精度磁强计、电场探测仪、GNSS掩星接收机、等离子体分析仪、高能粒子探测器、朗缪尔探针和三频信标发射机等8种载荷(Huang et al., 2018a).这里简单介绍本研究所涉及的用于电磁场观测的3种科学载荷的基本情况：

(1) 高精度磁强计(High Precision Magnetometer, HPM)由两个磁通门探头和绝对磁场校准装置组成，获取地球磁场矢量数据并实现对磁场总强度的测量，测量频率范围为DC~15 Hz(Zhou et al., 2018; Cheng et al., 2018; Pollinger et al., 2018)；(2)感应式磁力仪(Search Coil Magnetometer, SCM)提供10 Hz~20 kHz频段范围内的磁场波形和频谱数据，该载荷工作模式有详查、巡查、定标和数传测试四种工作模式，探头安装在一根4.5m长伸杆的末端，提供ULF(10 Hz~200 Hz)、ELF(200 Hz~2.2 kHz)和VLF(12.5 kHz~20 kHz)三个频段的变化磁场观测数据，采样率分别为1024 Hz，10.24 kHz，51.2 kHz(Cao et al., 2018); (3)电场探测仪(Electric Field Detector, EFD)由4根4.5 m长的伸杆向卫星本体外伸出的4个传感器组成，通过4个传感器之间的电位差换算出空间三维电场，提供ULF(DC~16 Hz)、ELF(6 Hz~2.2 kHz)、VLF(1.8 kHz~20 kHz)和HF(18 kHz~3.5 MHz)四个频段的电场数据，这四个频段的采样率分别为128 Hz，5 kHz，50 kHz和10 MHz(Huang et al., 2018b).

顶部电离层400~1000 km空间内另一颗重要的电磁卫星是DEMETER(Detection of Electro-Magnetic Emissions Transmitted from Earthquake Regions)，它于2004年6月发射入轨，2010年12月结束使命.DEMETER是世界上首颗科学目标专用于地震、火山及人类活动引起的电离层扰动研究的卫星(Parrot et al., 2006b).DEMETER卫星一直连续观测了6年多，在660~710 km的轨道空间内积累了翔实的观测数据.DEMETER没有搭载磁通门磁力计，因此无背景磁场观测数据.DEMETER主要搭载了感应式磁力仪(Instrument Magnetic Search Coil, IMSC)，用来观测19.5 Hz~20 kHz的变化磁场(Parrot et al., 2006b)和电场探测仪(Instrument Champ Electrique, ICE)来观测DC~3.175 MHz的电场(Berthelier et al., 2006).两颗卫星在电磁场观测方面的性能指标如表 1.

本研究主要使用了CSES卫星的HPM观测的总磁场三分量矢量数据和变化电磁场观测设备SCM和EFD在ULF/ELF/VLF频段记录的频谱数据和波形数据，以及DEMETER卫星的ICE和IMSC提供的电磁场19.5 Hz ~20 kHz频段的单分量频谱数据和0 Hz ~1.25 kHz频段的波形数据.关于六分量的电磁场波形数据，DEMETER卫星只提供地震危险区或者实验装置上空短暂时间内，1.25 kHz以下的六分量波形数据，通常只维持几分钟时长.相对优势的是，CSES卫星提供了ULF/ELF频段全轨道波形数据，在地震危险区上空提供VLF频段的短时波形数据，这为开展空间电磁波传播特征分析提供了丰富的数据.

2 张衡一号记录的波形数据频谱分析 2.1 波形数据频谱分析算法

本研究主要对CSES卫星观测的电磁场ULF/ELF/VLF频段的波形数据采用快速傅里叶变换(FFT)，将卫星记录的波形数据从时域转换到频域.FFT算法的优势是利用离散傅里叶变换(DFT)的奇偶虚实等特性，把数据序列分解成一系列的短序列，通过合理组合等方式避免重复计算、减少乘法运算次数，求出这些序列的DFT，从而达到快速计算的目的.DFT定义式可表示为

(1)

其中

直接计算式(1)需要N²次复数乘和复数加运算，N很大时计算量将会非常大，我们可以利用快速傅里叶变换(FFT)进行高效计算.假设N是2的整次幂，即N=2^m，若不满足上述条件，可对已知序列补零{x(1), …, x(N), 0, 0, …}来增加序列长度，使其长度为L(L为2的整数次幂).首先将(1)的求和过程分为两部分：

(2)

可以看出

(3)

则(2)可以写成

当k = 2p=0，2，…，时

(4)

当k=2p+1=1，3，…，时

(5)

式(4)和(5)是基-2FFT方法的基本过程，表示将长度为N的序列的DFT分解为两个长度为N/2序列的DFT，计算两大约N/2次复数乘和复数加.由于我们把N补零到2的整次幂，计算效率显著提高.

功率谱计算的目的是根据有限数据，给出信号的频率成分分布描述，是在频域内提取淹没在噪声中的有用信息的分析方法.对于有限长信号功率谱，可以运用直接法(周期图法)并且利用FFT变换来计算，对于有限时间序列x₀, x₁, x₂, …, x_N-1，采样频率为f=1/Δt，功率谱计算公式为

(6)

直接法得到的功率谱估值精度较差，并且方差不会随数据量增加而减小，也就是说大量的数据依然不能改善精度，针对这些存在的问题，可以通过多种方法进行改进，我们应用的是Welch法，也是最常用的PSD估计方法之一.

(1) 将信号序列X(n)进行分割成若干段，允许数据有重叠：

(7)

(2) 计算每段数据加窗的周期图：

(8)

其中P表示时间窗w_i的“功率”：

(9)

(3) 对加窗周期图进行平均，即可得到Welch谱估计：

(10)

为了降低被估计的PSD值的方差，Welch方法在数据分段时允许数据之间有重叠，这样就增加了式(1)中被平均的周期图数，并且在计算周期图时引入了时间窗可以有效提高PSD值的分辨率.

图 1中前三幅子图显示了CSES卫星在2018年9月16日记录的磁场ELF频段的原始三分量波形矢量.图 1最后一幅子图是计算的磁场的功率谱密度，从中可以观察到600 Hz附近存在一个明显的截止现象，截止频率通常与质子回旋频率非常接近并且位于质子回旋频率以下，图中黑色曲线表示当地质子回旋频率.在南极极区截止频率以上可以观察到很强的电磁辐射，并且演化为截止频率以下的嘶声波.低纬度地区存在可能由闪电引起的哨声波.我们在下文会讨论质子回旋频率的计算方法并验证其可靠性.

2.2 质子回旋频率附近电磁波动

我们利用质子回旋频率附近的电磁波动来验证算法，因为通常质子回旋频率附近会出现明显的电磁波动，并且会趋向于往质子回旋频率变化的方向传播，因此通过质子回旋频率的记录情况可以直接验证算法的正确性.通过单粒子轨道理论，利用带电粒子的洛伦兹力公式可得到质子回旋频率：

(11)

(12)

将式(11)带入式(12)可得

(13)

f_pe为质子回旋频率，q为质子电荷量，v为速度，r为拉莫尔半径，m_p为质子质量，取m_p=1.6726×10^-27kg，q=1.6022×10^-19C.

在本次计算中，我们利用CSES卫星的高精度磁强计记录的总磁场强度代入B，进而计算出沿轨道空间的质子回旋频率大小，如图 2a所示.图 2a显示了CSES卫星在2018年10月1日01:35—02:10 UT时间内记录的ELF 200~800 Hz的电磁波动情况，可以看出在600~700 Hz质子回旋频率附近，电磁波动主要沿着质子回旋频率附近出现强烈辐射.我们选择了DEMETER卫星日侧下行轨道的观测数据，通过对比分析可知CSES卫星和DEMETER卫星对相同频段的电磁场观测具有良好的一致性，结果如图 2b所示

从图 2可以看出质子回旋频率随着纬度减小而逐渐减小，这是因为地磁总场B在南北两极强，赤道附近较弱的原因造成的.CSES卫星的轨道空间在507 km左右，DEMETER卫星2010年轨道高度为660 km左右.对比两颗卫星分析结果显示：两颗卫星显示的质子回旋频率曲线符合基本规律，即两极高赤道低；在临近轨道空间分别在507 km和660 km高度上两条质子回旋频率曲线整体形态趋于一致；两颗卫星都在质子回旋频率附近观测到了明显的电磁波动，并且电磁波动在质子回旋频率附近有明显的截止现象.这也证明了CSES卫星与DEMETER卫星观测具有良好的一致性，也验证了我们的计算结果是可靠的.

2.3 闪电引起的ELF哨声波频谱特征分析

张衡一号卫星轨道空间存在着源自大气层的闪电或者雷暴引起的强烈电磁辐射，研究表明闪电产生的哨声波普遍存在于等离子体层内，这种电磁扰动主要表现为ELF/VLF频段的哨声波模.前人研究认为闪电产生的ELF/VLF频段范围内的脉冲辐射沿着地球—电离层波导管传播，并且能够穿过电离层(曹晋滨等，2009)，被近地轨道卫星观测到，通常卫星观测到的闪电事件表现为从几Hz到几十kHz瞬时增强，或者表现为高频信号先到，低频信号后到的降调哨声波模(Santolík et al., 2009; Inan et al., 2007).

DEMETER卫星曾观测到大量与闪电活动有关的哨声波(Santolík et al., 2009; Chum et al., 2009)，通常表现为降调模式的哨声波模，我们在CSES卫星也观测到了大量的类似降调模式的事件，但是更有趣的是发现了如图 3所示的上声调的哨声波事件，图 3是CSES卫星2018年10月4日15:53—15:56 UT时间内在赤道上空记录的ELF 200~400 Hz频段的哨声波动现象.

3 张衡一号记录的电磁波矢量分析 3.1 奇异值分解算法

关于电磁波的波矢量分析方法主要有Means、McPherron方法、Samson-Olson方法以及奇异值分解法Singular Value Decomposition(SVD)，其中SVD方法利用了谱矩阵的全部信息，对谱矩阵进行奇异值分解，获得波矢的最小二乘解，因而可以得到波场更加丰富的特性.本文主要采用的是Santolík等(2003)基于电磁卫星观测开发的电磁场奇异值分解(SVD)法，并将算法实现于张衡一号卫星电磁波动分析中.下面对该算法进行详细介绍：

若假定电磁波是单频平面波，则由Maxwell方程可知，波矢、电场矢量和磁场矢量两两垂直.对于给定的频率，有多种可行的频谱分析方法可以用来估算磁场和电场，例如，我们可以利用由快速傅里叶变换算法(FFT)实现的离散傅里叶变换(DFT)，并利用不同的加权窗口对实际样本进行时间序列分析，得到磁场的复数矩阵，利用可以得到一个Hermitian功率谱矩阵:

(14)

下标i和j表示磁场的三个笛卡尔分量，*表示共轭函数，〈〉表示在连续时间区间上的均值，其中是从实际样本的有限时间序列计算得到.由麦克斯韦方程磁通连续性定理得到磁场的散度为零:

(15)

如果方程(15)分别乘以B^*的三个笛卡尔分量，就会得到三个相互不独立的方程：

(16)

定义B_iB_j^*为信号功率谱矩阵S_ij：

(17)

则(15)中的三个复数方程可写成

(18)

其中包含了频谱矩阵S的所有分量以及波矢量列向量k.我们将方程重新整理成一个由六个方程组成的齐次方程组：

(19)

A是一个由谱矩阵的实部和虚部叠加组成的实矩阵：

(20)

由于(19)式右侧为零，两侧同时乘以任意实数不会改变等式结果，因此我们利用磁场只能计算波矢量的方向而无法求出波长，并且波矢量k乘以一个负数等式同样成立，但是所求波矢量k方向与矢量实际方向相反.因此我们只能在半个球面上定义波矢量方向.波矢量参数在球面坐标系中定义为

(21)

其中θ为极角(0°~90°)，ϕ为方位角(-180°~180°)，κ₁, κ₂, κ₃是归一化波矢量κ=k/k的三个笛卡尔分量.因为κ全部由θ和ϕ确定，其模始终为1.

方程组(19)中方程数量大于未知数数量，是超定方程组，从中任选两个独立的方程就可确定θ和ϕ的唯一解，Means方法(Means, 1972)、Sasmson-Olson方法(Samson and Olson, 1980)等谱分析都是基于此方法，尽管在表达形式上可能有所不同.

对于理想功率谱矩阵分析，利用(19)中的部分方程进行计算是可行的.而在实际应用中对于一个有限的时间序列，我们无法通过谱分析而完全确定功率谱矩阵，并且波通常包含未极化成分与噪声，这种情况下我们采用这种降级处理往往会丢失功率谱矩阵中的部分信息.SVD方法正是为了解决这个问题，用替代A，与之相对应的替换S，方程(19)为

(22)

如果方程没有简并为三个方程以内，就不可能得到同时满足这六个方程的.而利用最小二乘解可以找到一个使模最小的解，这就可以在不必实际解决最小化问题的情况下，只采用矩阵的SVD分析.

(23)

其中U是具有标准正交列向量的6×3矩阵，W是由3个非负特征值组成的3×3对角矩阵，V^T是具有标准正交行的3×3矩阵.

令w₁，w₂，w₃表示三个奇异值，并且满足w₁ < w₂ < w₃的关系，定义波磁场的平面度：

(24)

F趋近1，表示有1个或2个非零特征值，式(22)可以简化为一个或两个方程.有2个非零特征值时，可以用上述SVD方法得到唯一确定解.F小于1时，表示有3个非零特征值时，可通过求取最小二乘解来确定θ和ϕ.F值越接近0表示极化程度越低，F值为1是代表完全极化.定义极化椭圆的两个极化轴为

(25)

L_p等于或接近1，表示极化种类为圆极化，L_p等于或接近0，表示极化种类为线极化，并且只有1个非零特征值时，方程(22)可以简化为一个等式，未知数大于方程数，不能得到θ和ϕ的唯一解.

以上基于磁场三分量的奇异值分解，还有基于电磁场六分量的奇异值分解，具体详见Santolík等(2003)，在本文中略去.

3.2 电磁波的坡印廷能流算法

由于如上波矢量分析算法的局限性，我们无法确定波矢量k的方向，即无法分辨两个反向平行的波，为此我们利用坡印廷(Poynting)矢量的方向和大小来辅助解决这一问题.坡印廷矢量即能流密度矢量，方向为电磁能的传播方向，大小表征单位面积的能量传输速率(J·m^-2·s^-1)，若电场强度为E，磁场强度为H，则该处能流密度定义为

(26)

S即坡印廷矢量，传播方向沿电磁波传播方向，由E和H按右手螺旋定则确定.电、磁场强度的变化引起能流密度的变化，通过对变化电磁场中坡印廷矢量方向的研究就可以确定波动的传播方向.利用六分量电磁场数据我们已经求取了6×6的谱矩阵，利用谱矩阵可以通过下式计算平均坡印廷矢量：

(27)

其中E(f)和B(f)分别表示电、磁场波形的傅立叶变换，μ₀为真空磁导率，将上式各分量展开为

(28)

其中为电磁场互功率谱的实部，S_Sk是Poynting矢量(i, j, k=1, 2, 3)笛卡尔分量上有效值的谱密度.

3.3 场向坐标系统构建

我们定义电磁波波矢量为k，并且k垂直于电磁波的极化平面，根据麦克斯韦方程

(29)

可知由电场矢量E和磁场矢量B可以得到波矢量k.其中ω=2πf为电磁角频率.定义坐标系如图 4所示，B₀为地磁总场方向，X轴在磁子午面内由地心垂直于地磁线指向地球外部，Y轴垂直于磁子午面，方向由右手定则确定，Z轴即地磁总场B₀的方向，波矢量k与地磁场B₀的夹角为θ，ϕ为波矢量k在磁子午面上的投影与X轴的夹角.由方程组(29)可知磁场矢量B总是垂直于波矢量k和电场矢量E，这一性质使我们有时域和频域两种方法来确定波矢量.时域分析常用的方法为磁场扰动的最小变量分析法，频域分析为了更真实的反映波的信息通常采用电磁场多分量信息进行分析.

首先将CSES的高精度磁强计HPM记录的总磁场三分量构建B₀数组，将感应式磁力计SCM记录的磁场三分量构建B_w数组，电场仪EFD记录的电场三分量构建E_w数组.并将这三组数据从卫星球坐标转换成直角坐标系.通过B₀构建标准化的场向坐标系，其Z分量代表总磁场的强度.然后将B_w和E_w转换到这个场向坐标系中，根据3.1节和3.2节的算法逐步计算出波矢量的传播参数.其具体结果将结合第4节中的典型波动事件进行详细描述.我们已将本算法应用于DEMETER卫星观测的波形数据分析中，验证了算法的正确性，在此不再显示其对比结果，具体计算结果下一节将结合典型电磁波动事件介绍.

4 张衡一号卫星记录的典型电磁波动事件的波矢量分析

张衡一号卫星自发射以来，在电离层顶部507 km左右的空间记录了丰富的电磁波动.这些波动的频率覆盖了从几Hz至几个MHz的频率范围.在此我们挑选了两个典型波动事件，利用前文介绍的算法开展了电磁波动频谱分析及波矢量计算.

4.1 ELF/VLF准周期辐射波

图 5是CSES卫星在2018年9月17日08:06—08:11 UT时间内记录的ELF/VLF频段内一个准周期辐射现象.其中图 5a显示了感应式磁力计记录的变化磁场三分量波形转换的总功率谱密度(Power Spectral Density)，可以看出卫星记录到了很明显的离散结构的波谱结构.这些波谱之间的周期不固定，如图 5b显示了其在ELF频段1200 Hz的功率谱密度变化情况，其周期变化范围从大约16 s至30 s不等.因此学术界将其定义为准周期辐射(Quasi-Periodic emissions)(Sato et al., 1974; Hayosh et al., 2016).根据前人对DEMETER卫星记录的准周期辐射的研究(Hayosh et al., 2016；Němec et al., 2016)，可以发现CSES卫星记录的准周期辐射具有清晰的离散结构.

准周期辐射是一种重要的电磁波动现象，其产生机制比较复杂，一种产生机制与ULF地磁脉动有关，另一种产生机制与能量粒子相互作用时造成的回旋不稳定性有关(Hayosh et al., 2013；Němec et al., 2016).本文主要关注点为张衡一号卫星电磁波动分析方法的实现和验证，因此在此不再深入讨论其产生机制.我们继续利用SVD波矢量分析算法，计算出了该波动事件的传播参数，如图 6所示.

图 6显示了2018年9月17日08:04—08:11 UT时间范围内利用SVD方法计算出的波传播特征.图 6a和6b分别显示了变化磁场、电场的功率谱密度变化情况，可以看出同时在电场、磁场显示出了准周期辐射现象，证明这个波动是电磁波动.图 6c显示椭圆率(Ellipticity)的变化情况.根据上文3.1节SVD算法可知，椭圆率接近1表明该波是右旋圆极化，-1表示左旋圆极化，0表示线性极化.可以看出本波动事件的椭圆率接近1，表明该波动是右旋极化的哨声模波动.图 6d表示了波矢量k与总磁场B₀的夹角传播角θ大概在30 °左右变化，说明波是沿着磁力线斜向传播.图 6e显示了波矢量k的方位角ϕ大概为180°，证明波是往地球方向传播，说明其来源于高于卫星轨道空间的地方而非地球方向.图 6f表示波动的平面度，平面度越高表示利用SVD方法所得结果的可信度越高.图 6g和h分别显示了Poynting能流与B₀的两个夹角，即极角θ_s和方位角ϕ_s，可以发现Poynting能流以平行于B₀方向为主.

Hayosh等(2016)根据DEMETER卫星的观测对该类型的电磁波动进行了SVD计算并对将近200个同类事件的波矢量参数进行了统计分析，获得的结果认为这种类型的右旋电磁波动通常沿地磁场斜向方向传播，其波矢量变化范围为30°~70°，通过与之比较，本计算结果符合准周期辐射波动的典型特征.

4.2 ELF/VLF嘶声波

电离层嘶声波(hiss)是电离层中频率范围为100 Hz至3 kHz的一种普遍存在的宽频带、无结构的电磁波(Chen et al., 2017; Zhima et al., 2017)，通常表现为右旋极化模式，且嘶声波表现出了无结构的波谱特性，有研究者认为这种无结构嘶声波可能来源于结构明显的合声波(Santolík et al., 2006; Parrot et al., 2016)，合声波在传播过程中结构逐渐散失，从而表现出嘶声波无结构的频谱特征(Zhima et al., 2013).前人研究发现即使在地磁平静日嘶声波仍然广泛存在，但在太阳活动引发的磁暴、亚暴等事件发生期间嘶声模辐射会加剧(Smith et al., 1974; Thorne et al., 1974; Meredith et al., 2004).

Zhima等(2017)通过研究DEMETER和THEMIS卫星的共轭观测，发现电离层嘶声的产生源可能有三种：(1)较高轨道空间(等离子体层、内磁层)产生的等离子体嘶声下传到电离层高度，形成电离层嘶声并被低轨道卫(DEMETER，CSES)记录到；(2)大气层产生的闪电事件也可能是电离层嘶声的产生源(embryonic source)，通常闪电触发的嘶声表现在2 kHz以下；(3)低轨道卫星空间本身等离子体不稳定性也可能产生嘶声波动.然而经过近半世纪的卫星观测研究，电离层嘶声的产生机制依然尚未真正揭开，仍旧是学界研究热点.

图 7是CSES卫星2018年10月1日在轨道号0036631上记录的ELF频段辐射的嘶声现象，嘶声事件发生在05:32—05:41 UT位于质子回旋频率以下200~400 Hz处，表现为无结构状的密集电磁辐射.

我们提取了05:32—05:36UT时间范围内的6分量电磁波形数据，利用奇异值分解(SVD)法进行波矢量分析，其结果如图 8所示.图 8a—b显示了变化磁场、电场中波的频谱特征，图 8c显示了椭圆率，其值大概趋近于-1，说明这次波动为左旋极化，这很大可能是由离子回旋运动产生，经过计算氦离子在这个轨道空间的回旋频率大概在200 Hz左右，因此这个左旋的嘶声极化特征是可能存在的.图 8d显示了波矢量k传播极角θ约为30°，图 8e显示了波矢量k的方位角为70°~80°，说明波斜向传播.图 8f中平面度并不是全部接近1，说明这个嘶声波传播的方向并不是集中在一个方位上.

5 讨论及结论

本研究主要利用张衡一号电磁卫星记录的ULF/ELF/VLF的电磁场波形数据，研发了电磁波波矢量分析工具，该工具包括了电磁波形频谱变换和基于奇异值分解(SVD)方法的波矢量分析算法.基于奇异值分解SVD算法考虑了谱矩阵中包含的所有信息，能够确定电磁波的振幅、相位、传播角、频移、极化等参数.SVD方法只给出了半球中波矢量±180°的方位，为了解决±180°的方位偏差，我们在工具中增加了Poynting矢量分析算法，辅助SVD算法确定波动的波矢量方向.由于SVD算法会计算每一个采样点，因此如果观测数据较长，将消耗大量的计算时间，通常全轨道数据可能花费几个小时的CPU时间.因此建议在进行SVD计算前，先单纯绘制频谱图，根据频谱图大概确定电磁波动的起始和结束时间，然后进行针对波形事件范围内的SVD波矢量分析，不建议全轨道直接计算.其次时间校准问题，通常情况下CSES的载荷开机时间不是严格同步，因此需要在数据预处理前把HPM，EFD和SCM的数据根据最晚观测时间进行截取，使得计算的波形事件具有同步的电磁场观测数据.最后，由于每个典型波动事件的强度不同，背景电磁噪声也不同，因此在计算结果绘制过程中，需要设置一定范围以下不显示，才能突出典型事件的波矢传播特征.

关于波矢量分析方法典型代表有Means，McPherron方法，Samson-Olson方法、奇异值分解法Singular Value Decomposition(SVD).相对来说，SVD方法利用了谱矩阵的全部信息，对谱矩阵进行奇异值分解，获得波矢的最小二乘解，因此能够得到更加丰富的波场信息.在研制本软件包过程中，我们也采用了Means方法进行校验结果的准确性.在研究电磁波的传播特性中，了解能流的方向和大小非常重要.通常除了波矢量分析获得波的传播角以及极化参数以外，我们还需要了解波的能流，即坡印廷(Poynting)矢量的方向及大小.由于采用单一磁场SVD方法分析波矢量时，存在一定的不确定性，即两个反平行的波矢量方向的确定，因此需要坡印廷矢量分析来辅助判断波矢方向.SVD算法是用于电磁多分量数据处理的一种重要和可靠的方法，但是此方法只对平面波有意义.如果波场性质更加复杂，例如远处同时存在两个辐射源，不满足单频平面波的假定，则波矢分析必须采用波分布函数(Wave Distribution Function，WDF)方法.因此未来我们还将考虑WDF算法的实现，使此软件包更加完善.

本文通过对张衡卫星观测的典型电磁波动事件如闪电引起的哨声波、准周期辐射、电离层嘶声波等进行频谱分析和矢量计算，验证了波矢量分析工具的正确性，并通过相同事件相同算法对DEMETER卫星观测结果的对比分析验证了程序的正确性.通过与DEMETER卫星的对比研究，也证实了张衡一号电磁卫星的观测与DEMETER卫星具有良好的一致性，而且张衡一号卫星的电磁场采样频率相比于DEMETER卫星更高，更能记录到精细的波谱结构.如本次研究中提及的准周期辐射波动，张衡一号卫星记录了更清晰的波包结构，这方面我们进行了深入分析，相关结果正在撰写新的研究论文.

顶部电离层空间是一个各圈层紧密联系的复杂耦合系统，直接受太阳活动、行星际扰动、岩石圈，大气层，内磁层的影响.从上而来的扰动主要源于太阳的强烈活动，如太阳耀斑、日冕物质抛射、共转相互作用区等活动引起的近地空间等离子体、电磁场剧烈扰动(磁暴/亚暴).从地球方向而来的扰动，主要包括地震、闪电、台风、大气潮汐、人工甚低频发射站等引起的各种尺度的电磁波动.张衡一号电磁卫星将在顶部电离层空间长时间持续观测，预期在轨运行可至2025年，为研究这个复杂的耦合系统带来重要的第一手观测资料.通过前人的研究和本研究经验，可知这个空间内的电磁波的产生机制、传播特征以及来源也各有不同，同时太阳活动和行星际扰动也加剧了卫星轨道空间电磁波动的复杂性.因此利用电磁卫星数据进行电磁波动分析时，应考虑多种因素的干扰，更准确的将电磁波信号与干扰信号区分开来.未来我们将持续利用该方法对张衡一号卫星记录的电磁波动事件开展深入的研究，探索其传播和产生机制.

致谢  本工作使用了中国国家航天局和中国地震局支持的张衡一号卫星的观测数据(http://leos.ac.cn)和法国DEMETER卫星观测数据(https://demeter.cnes.fr/).本研究由国家自然基金面上项目(41874174，41574139)和亚太空间合作组织地震研究项目(APSCO Earthquake Research Project Phase II)资助.感谢北京航空航天大学符慧山教授对SVD算法的帮助和指导.