2021, Vol. 64
2. 应急管理部国家自然灾害防治研究院, 北京 100085
2. National Institute of Natural Hazards, Ministry of Emergency Management of China, Beijing 100085, China
利用电磁卫星监测地震始于20世纪80年代，现已探测到了大量的电磁异常信息(Shen et al., 2018; Yan et al., 2018; Larkina, 1983; Chmyrev et al., 1989; Parrot and Mogilevsky, 1989; Pulinets and Legen′ka, 2003; Molchanov, 1993; 蔡军涛等, 2007)，从而为地震电磁扰动提供了新的观测手段.然而，引起电磁异常扰动的因素很多(尼鲁帕尔·买买吐孙和张永仙, 2012)，大量研究结果显示，低层大气气象活动中的一些剧烈的天气现象也会引起电离层扰动，譬如台风、飓风、寒潮、龙卷风、闪电等(刘依谋等, 2006; 肖赛冠等, 2007).在众多的气象活动中，闪电是发生最频繁的一种自然现象，Christian等(2003)使用光学瞬态探测器探测到全球闪电频数为44 fl/s(fl为flash简写，表征闪电发生的次数).闪电放电过程中激发宽频带的电磁波，频率范围在500 Hz~10 kHz，因类似于口哨的声音，故被称为哨声波(Helliwell, 1965).
闪电哨声波能够到达电离层、磁层，是探测空间物理环境的重要媒介(Helliwell, 1965).Záhlava等(2018)通过分析DEMETER卫星和RBSP卫星的闪电哨声波发现等离子体层内部的哨声波模的纵向依赖性很强.Bayupati等(2012)通过分析AKEBONO卫星的闪电哨声波的色散发现闪电哨声波的色散趋势是确定等离子体层中整体电子密度分布的有力方法.Oike等(2014)分析AKEBONO卫星探测的闪电哨声波发生频率证明了电离层中的闪电哨声波与闪电活动、地球周围的电子密度分布密切相关.Clilverd等(2002)借助全球闪电定位网(WWLLN)来确定闪电的来源位置，再根据闪电哨声波从来源点到观测点的传播时间，最后利用传播理论导出电子密度沿传播路径的函数.可见了解闪电哨声波的物理参数和位置等是开发利用闪电的基础，是研究电离层和磁层等空间环境的重要技术手段.
闪电哨声波观测主要有地基观测和天基观测.20世纪80年代，我国开始探索地基闪电哨声波观测，利用自主研发的地面哨声波接收机对低纬哨声波传播特性进行大量的研究(王挺柱等, 1982；徐继生等, 1989；夏利东和王友善, 1994；田茂等, 1991)，但受限于当时模拟电子技术与计算机技术，哨声波接收机灵敏度不高，不易存储(Chen, 2017).武汉大学于2016年成功研制新型甚低频接收机系统，已在我国低纬地区进行多站点(随州、乐山、武汉)联合观测，取得了突破性进展(Chen et al., 2016, 2017; Zhou et al., 2020; Yi et al., 2020).
2018年2月，我国首颗电磁卫星张衡一号(ZH-1)发射成功，具备了天基观测闪电哨声波的能力.ZH-1卫星覆盖南北纬65°，在中国大陆及周边1000 km区域及全球两个地震带(太平洋地震带和欧亚地震带)进行详查模式的观测，其他区域为巡查模式.ZH-1卫星在轨飞行高度约507 km，其位置接近电离层顶部和等离子层边界，这个区域有丰富的ELF/VLF频段波动事件，如闪电哨声波、准周期辐射等(Zhima et al., 2020).ZH-1轨道倾角97.4°，属于太阳同步轨道，降交点地方时为下午2∶00；轨道回归周期为5 d，即每5 d星下点轨迹相同；在一个回归周期内能够实现全球约500 km空间分辨率的观测(袁仕耿等，2018).卫星绕地球飞行一圈约94 min，大部分载荷在±65°纬度范围开机工作，观测数据按升轨(夜晚)和降轨(白天)分别存储，每半轨(升/降轨)观测约34 min；在同一天内相邻的升轨(或降轨)空间分辨率约2000 km.所搭载的感应磁力仪载荷(SCM)，通过法拉第电磁感应定律获得电离层的感应磁场数据，能够捕获全球闪电哨声波信号，其在巡查模式下仅获得功率谱数据.到目前为止，ZH-1已经在轨观测3年多，采集了大量的全球电磁场的波形和功率谱数据，其中SCM的3分量X/Y/Z包含3个频段ULF/ELF/VLF，频点范围ULF: 10~200 Hz，ELF：200 Hz~2.2 kHz，VLF：12.5 Hz~20 kHz，原始数据的采样率为51.2 kHz，功率谱数据的频点间隔ULF: 0.25 Hz，ELF: 2.5 Hz，VLF：12.5 Hz，详查模式VLF波形数据80 ms包含4096个点(Wang et al., 2018; Zhou et al., 2018)，每天产生大约10 G的数据量.如何从如此庞大的观测数据中自动识别闪电哨声波事件尤为关键和紧迫.
目前基于海量闪电哨声波数据开展全球时空分布规律和相关参数的研究较少，主要存在以下两个方面的原因：一是需要足够的闪电哨声波事件；二是缺少有效的闪电哨声波自动识别算法，依靠人工从大量的电磁观测数据中获取闪电哨声波事件是一件极具有挑战性的工作.比如：Zhou等(2020)针对武汉VLF地面观测数据中隐藏的闪电tweek现象，通过设置能量谱阈值和时间宽度阈值的方式提出了简单快捷的自动识别算法.Yi等(2020)借用该算法对随州站2016年2月3日至2月29日的WHU ELF/VLF接收机数据中的所有tweek现象进行自动识别和统计分析，发现该站tweek事件的日发生率的差异很大，数量大约在800至6000个之间.然而，其采用的识别算法(Zhou et al., 2020)的误报率较高，主要原因是基于阈值匹配的分类策略使干扰噪声对识别结果的影响较大，该方法明显不适用于电磁卫星数据.电磁卫星的观测数据中携带了复杂多样的空间电磁等扰动信号，如闪电哨声波、嘶声、VLF地面发射源、磁暴等现象.现有的研究主要是集中在一些明显的太阳活动、地磁活动及VLF地面发射源对电离层的扰动方面(李柳元等2011; Li et al., 2015; Yang et al., 2020; 廖力等，2019; Zhao et al., 2019)，较少有研究如何准确地从观测数据中自动识别闪电哨声波.随着电磁观测大数据池的形成，自动识别闪电哨声波是自动提取闪电哨声波轨迹信息及其相关物理参数的重要前提，是对空间物理环境进行全面深入分析的关键环节.闪电哨声波源于闪电激发的宽频带电磁波，其高频和低频成分的速度差使其在频谱图上呈现明显的“色散”状特征.当路径长度较长或沿传输路径的电子密度较密集时，色散变大，可见色散是反映等离子体层中的电子密度的重要参数.Helliwell(1965)根据地面VLF观测数据按照其色散特征将闪电哨声波的形态分为9类，详细内容如表 1所示，值得注意的是“多成分”类型中存在两种不同的子类型.Santolík等(2009)指出地面VLF观测站记录的闪电哨声波的形态与卫星观测到的形态极其相似且均具有频散特性，如图 1和图 2所示.
|
表 1 闪电哨声波分类(Helliwell, 1965) Table 1 Type and principal definitions of lightning whistlers (Helliwell, 1965) |
|
图 1 AKEBONO卫星闪电形态图例(Bayupati et al., 2012) Fig. 1 Lightning whistlers of AKEBONO (Bayupati et al., 2012) |
|
图 2 DEMETER卫星闪电形态图例(Parrot et al., 2015, 2019) Fig. 2 Lightning whistlers of DEMETER (Parrot et al., 2015, 2019) |
目前基于电磁卫星闪电哨声波自动识别的研究主要是针对AKEBONO卫星采集的波形数据(Bayupati et al., 2012)、Arase卫星采集的数据(Ali Ahmad et al., 2019)以及DEMETER卫星采集的频谱数据(Parrot et al., 2019).
按照Eckersley′s的定理(Gurnett et al., 1990)，闪电哨声波的时间和频率之间满足公式(1):
|
(1) |
其中D是闪电哨声波色散常量; t是以频率f到达的时间，t0是闪电触发的时间.Bayupati等(2012)将该定理应用到AKEBONO卫星的闪电哨声波识别中.传统的闪电哨声波自动识别方法是通过识别功率谱中的直线的方式实现的，如图 1所示：t1是在频率f1到达的时间，t2是频率f2到达的时间，D的计算公式如式(2)所示：
|
(2) |
图 1显示闪电哨声波是有一定斜率的一条直线，而公式(1)(2)仅仅表示了闪电哨声波频率和时间之间的简单关系.
Ali Ahmad等(2019)认为上述方法很难检测到Helliwell(1965)总结的其他闪电哨声波类型，重新整理了Arase卫星的闪电哨声波的类型，并根据不同的类型设计了不同的模式特征.其识别的方法是：首先用短时傅里叶变换(STFT)计算频谱，接着采用高斯卷积核对频谱图进行模糊处理，再用拉普拉斯算子做边缘提取，然后用图像分割算法分割图像，如图 3所示；并将图像平均划分为若干个格子，如图 4所示；最后，采用Bresenham′s直线检测算法提取特征(Bresenham, 1977)并输入到决策树进行识别，识别精度只有75%.该算法存在以下几个问题：(1)识别过程中需要大量的图像预处理操作，比如高斯卷积核的参数选择、拉普拉斯算子参数的选择、分割算法的参数选择等，这些参数的选择对采用Bresenham算法提取特征的影响较大；(2)Bresenham算法中需要设置参数，该参数的选择会影响特征的鲁棒性；(3)该算法的评测指标只采用精度作为评价指标，该指标无法评估分类器的召回能力和排序能力.
|
图 3 分割结果(Ali Ahmad et al., 2019) Fig. 3 Segmentation results (Ali Ahmad et al., 2019) |
|
图 4 格线法(Ali Ahmad et al., 2019) Fig. 4 Grid (Ali Ahmad et al., 2019) |
总之，目前的电磁卫星闪电哨声波识别算法存在下列问题：
(1) 模式特征过于简单导致其鲁棒性差；
(2) 分类技术粗糙导致其泛化能力不强；
(3) 以AKEBONO卫星和Arase卫星的数据为主，基于ZH-1卫星的闪电哨声波(如图 5所示)的识别研究未见开展.
|
图 5 张衡一号卫星的SCM功率谱数据的时频图中闪电形态 Fig. 5 Lightning whistlers from ZH-1 satellite |
综上，本文在ZH-1卫星SCM数据的基础上，探索基于机器学习的闪电哨声波自动识别框架，根据闪电哨声波色散形态创建L形态卷积核提高特征鲁棒性，采用支持向量机技术提高分类的泛化性能.
1 数据收集数据主要来自张衡一号卫星在2019年8月份的三个轨道的SCM载荷的VLF波段的波形数据，根据资料：大部分的闪电哨声波在100~6000 Hz区域内存在L形态，且持续时不会超过20 s(Fiser et al., 2010)，因此，本项目的数据收集过程如图 6所示：设计20 s的时间滑动窗，步长是4 s，截取波形数据，做短时傅里叶变换得到时频图，并创建闪电哨声波数据集.该数据集含316个闪电哨声波时频图像，8000个非闪电哨声波时频图像.由于本研究的重点是识别时频图像中是否具有色散形态特征，因此时频图中未绘制横坐标、纵坐标和色阶栏等辅助信息，其样本示例如图 7所示.
|
图 6 数据采集 Fig. 6 Data acquisition |
|
图 7 样本示例 (a) 闪电哨声波样例；(b) 无闪电哨声波样例. Fig. 7 Sample example (a) Lightning whistlers; (b) No lightning whistlers. |
闪电哨声波识别算法的目标是识别每个时频图中是否有闪电哨声波.目前深度学习模型是机器学习的主流模型，但仍需要大量的数据样本作为基础.而本文目前收集的闪电哨声波数据量较少，基于此仍然用传统的识别算法技术探索闪电哨声波识别方案，方案包括训练过程和识别过程，如图 8所示：实线是训练流程，虚线是识别流程.训练流程的主要目的是获得闪电哨声波的识别模型，识别流程的主要目的是应用闪电哨声波识别模型.训练过程包括：灰度化处理、尺度处理，模糊卷积处理，L形态卷积处理，训练SVM模型.其识别流程包括：灰度化处理、尺度处理，模糊卷积处理，L形态卷积处理，应用SVM模型进行闪电识别.
|
图 8 识别算法(实线是模型训练流；虚线是模型识别流) Fig. 8 Recognition algorithm (the solid line represents training process; the dotted line denotes recognition process) |
观察原始的样本频谱图，如图 9所示，发现识别闪电哨声波的依据是频谱图中是否存在L形态特征，无需颜色特征，则首先采用灰度化处理消除颜色的影响.对图 9按照公式(3)进行灰度化处理(章毓晋, 2016)，得到如图 10所示的灰度图.
|
(3) |
|
图 9 彩色频谱图(有闪电) Fig. 9 Colorful spectrum (lightning whistlers) |
|
图 10 灰度化(有闪电) Fig. 10 Grayscale (lightning whistlers) |
其中RGB表示原始的样本频谱图，RGB.R是其红色通道的像素值，RGB.G和RGB.B分别是绿色和蓝色通道的值，Gray是灰度化的频谱图.
2.2 尺度处理图像的尺度越大凸显的细节越明显，但整体的形态特征不够突出，如图 11所示；而尺度越小，容易突出其整体的形态特征，本文将频谱图像尺度缩小到50×50，得到如图 12所示的图像.对比图 11和图 12：发现图 12的L形态特征较图 11更为明显，且降低的数据维度更有利于计算.
|
图 11 原尺度 Fig. 11 The original scale |
|
图 12 小尺度 Fig. 12 The smaller scale |
观察图 12发现该图像条纹状纹理比较严重，表明该图具有明显的阶跃边缘.这些边缘信息不利于L形态特征的识别.为了弱化阶跃边缘对识别的影响，本文设计如图 13a所示的模糊卷积核，用其对图 12进行卷积处理，得到如图 13b所示的效果：其阶跃边缘信息被弱化了.
|
图 13 模糊卷积 (a) 模糊卷积核；(b) 模糊卷积效果. Fig. 13 Fuzzy convolution (a) Kernel; (b) Result. |
为了进一步增强闪电哨声波L形状特征提高分类性能，根据L形态的特点，设计如图 14a所示的9×9的L形态卷积核(关于卷积核尺度为何是9×9，将在讨论部分进行详细分析)，采用该卷积核对图 13b进行卷积计算，得到如图 14b所示的效果，增强了L形态特征.
|
图 14 L特征卷积 (a) L特征卷积核；(b) 卷积效果. Fig. 14 L characteristic convolution (a) Kernel; (b) Result. |
将图 14b的图像输入到支持向量机(SVM)分类器进行训练和识别.SVM是基于统计理论的监督分类器(Sánchez and David, 2003).本文中的类别主要有两种：闪电哨声波和非闪电哨声波.首先，将时频图数据分为训练集和测试集.
训练集的图像数据被表示为：x1，x2，…，xn，其相应的类别标签是y1，y2，…，yn，yi∈[-1，1]，其中-1表示该样本中无闪电哨声波，1表示该样本中有闪电哨声波.
线性SVM分类器的主要目标是找到一个超平面w=(w1，w2，…，wn)，能够将“-1”类和“1”类更好地分开.该超平面的数学模型如公式(4)所示：
|
(4) |
其中x是在超平面上的样本点，b是偏置.
若存在使得分类间隔最大的超平面，其满足下列条件：闪电哨声波的样本正确地落在闪电哨声波的分类区，则有公式(5)；
|
(5) |
非闪电哨声波的样本正确落在非闪电哨声波的分类区，则有公式(6)；
|
(6) |
公式(5)和公式(6)合成一个公式(7)：
|
(7) |
SVM将在众多的分类超平面中，自动寻找使得特征空间间隔最大的超平面.其数学模型表示为公式(8)：
|
(8) |
训练模型的目的是找到满足式(8)的w.
由于大部分情况下，数据无法满足线性可分的要求，此时满足上述条件的超平面不存在，对于这种非线性可分的问题，SVM提供核函数k(xi，xj)代替公式(5)，将原数据映射到高维空间，使得其在高维空间线性可分(Sánchez and David, 2003).本文中的核函数采用的是多项式核函数如式(9)所示：
|
(9) |
由于闪电哨声波的样本与非闪电哨声波的样本不均衡，按照传统的机器学习的惯例：正负样本的比例不能超过1∶3.而本文中的闪电哨声波样本是316个，与之对应的非闪电哨声波的样本数据量不能超过1000个，实际上本文获得的非闪电哨声波的样本是8000个.为了测试提出的闪电哨声波识别算法是否有效，制定下列实验方案.
1) 数据集
数据集由316个闪电哨声波样本和1000个非闪电哨声波样本组成.1000个非闪电哨声波样本是从8000个非闪电哨声波数据中进行随机选择而得到.该过程得到316个闪电哨声波组成的正样本集和1000个非闪电哨声波组成的负样本集.
2) 训练集和测试集
从正样本集随机取80%的数据作为训练样本集，同样从负样本集随机取80%作为训练样本，该过程得到含有正样本和负样本的训练集.
将正样本集和负样本集中剩下的样本组合为测试集.
3)特征提取
针对训练集，用以下三种不同的特征提取方法提取时频图像特征：
(1) 原始灰度图像特征，用Gray表示.
(2) 对原始灰度图像进行模糊卷积处理后的特征，用Gray_Blur表示.
(3) 先对原始灰度图像进行模糊卷积处理，再采用本文提出的L形态卷积核进行处理得到的图像特征，用Gray_Blur_L表示.
4) 训练
用三种不同的特征分别训练SVM模型，得到三种不同的基于SVM的闪电哨声波识别模型.
5) 测试
针对测试集的数据，提取三种不同的特征，并分别输入三个不同的SVM识别模型进行识别.
6) 评估识别效果
对识别效果采用四种指标进行评估：精度(Precision)、召回率(Recall)、F1值和ROC面积(AUC).
7) 重复实验
重复上述实验过程10000次.
3.2 SVM参数选择本算法中的SVM来自Python的SVC库.选择的核函数是多项式核，如公式(9)所示.通过交叉验证得到的最佳参数d是13，其他均使用该库下的默认参数.
3.3 评价指标采用四种不同的评价指标评估每一次实验的效果，设P表示正类样本，N表示负类样本，接下来将阐述四种指标的定义.
首先，符号定义如表 2所示.
|
|
表 2 符号定义 Table 2 Symbol definition |
(1) 精度(Precision)的定义如式(10)所示：
|
(10) |
含义：在预测出的所有正例中(TP+FP)，真正例(TP)占的比例，通常该指标越大越好.
(2) 召回率(Recall)的定义如式(11)所示：
|
(11) |
含义：在数据集中所有的正例中(TP+FN)，真正例(TP)占的比例，通常该指标越大越好.
(3) F1值(F1-Score)的定义如式(12)所示：
|
(12) |
其含义：通常情况下，Precision高，Recall就低，Recall高，Precision就低.指标F1-score，综合考虑Precision和Recall的调和值.Recall不变时，Precision越大，1/Precision越小，从而F1越大.同理：Precision不变时，Recall越大，1/Recall越小，从而F1越大.该指标1越大越好.
(4) AUC-ROC：ROC(receiver operating characteristic curve, ROC)曲线下的面积，该值越高意味着分类的排序性能越佳.
3.4 评价策略在每一次实验的训练集和测试集上，均使用三种不同的特征提取方法提取图像特征，并以此训练得到三个不同的分类器，最后采用精度、召回率、F1和ROC-AUC评估每一个分类器的性能.由于每次的训练集和测试集不同，单次的四个评估指标难以充分评价本文提出的闪电哨声波识别算法的效果，因此，开展实验10000次，并在四个评估指标的基础上制定了如下的评价策略：
(1) 部分识别结果展示.
(2) 总体识别效果评价：对10000次实验的四个评价指标值分别求平均值，以此来评价不同算法的总体识别效果.
(3) 识别效果稳定性评价：对10000次实验的四个评价指标采用盒形图评估识别算法的稳定性.
(4) 识别效果差异性评价：为了评价不同的特征之间导致的分类效果是否具有明显差异，采用T检验进行差异性评价.阈值p=0.05，即小于0.05为具有明显差异，若大于0.05表明不具有明显差异.
3.4.1 部分识别结果展示部分识别结果如图 15和图 16所示：图 15是识别闪电哨声波的部分结果，(a)是正确识别的时频图，(b)误识别的时频图像.之所以将图 15b识别错误的原因是：其L形态特征(实线箭头所示)相比右边的“嘶吼”(虚线箭头所示)特征明显微弱很多，该特点在灰度图中更为明显，如图 15c所示，此时时频图的主要信息体现在嘶吼的特征上，则分类器出现了误判；图 16是识别非闪电哨声波识别的部分结果，(a)是正确识别的结果，(b)是误识别的结果.识别错误的原因是图像中出现的多重“嘶吼”造成了L形态特征的假象，如图(c)(d)(e)所示，分类器出现了误判.
|
图 15 闪电识别 (a) 正确识别；(b) 错误识别；(c) 灰度化图(b). Fig. 15 Lightning identification (a) Correct identification; (b) Error identification; (c) The gray image of the (b). |
|
图 16 非闪电识别 (a) 正确识别；(b) 错误识别；(c) 灰度化图(b)；(d) 模糊卷积图(c)；(e) L形态卷积核处理图(d). Fig. 16 The performance of recognizing the image without lightning whistles (a) Correct identification; (b) Error identification; (c) Grayscale map (b); (d) Processing the (c) with fuzzy convolution; (e) Processing figure (d) with the L morphological convolution kernel. |
本小节主要针对10000次的Precision、Recall、F1和AUC-ROC求平均值，度量不同图像特征下的闪电哨声波分类效果，定量结果如表 3和表 4所示.
|
|
表 3 10000次实验后的平均效果 Table 3 The performance obtained by 10000 experiments |
|
|
表 4 ROC面积AUC的均值 Table 4 The mean of the AUC of the ROC area |
(1) 基于Gray_Blur_L特征的闪电哨声波分类效果，在Precision、Recall和F1score三个指标上明显优于Gray_Blur特征和Gray特征.比如观察表 3中的精度指标，Gray_Blur_L特征使闪电哨声波识别精度达到0.945，优于Gray_Blur特征的0.915和Gray特征的0.732.
(2) 基于Gray_Blur_L特征的闪电哨声波分类效果，在AUC指标上明显优于Gray_Blur特征和Gray特征.比如，Gray_Blur_L特征使得分类器的排序能力达到0.989，优于Gray_Blur特征的0.975和Gray特征的0.882.
总之，本文通过创建L形态卷积核提取的Gray_Blur_L特征在闪电哨声波分类的各项性能方面表现最佳.
3.4.3 识别效果稳定性评价为了对分类的稳定性进行定性评价，根据每个指标的10000个数值绘制如图 17和图 18的盒形图.
|
图 17 精度、召回率和F1的盒形图 (a) 闪电哨声波；(b) 非闪电哨声波. Fig. 17 The box diagram of the precision, recall rate and F1 (a) Lightning-whistler; (b) Non-lightning-whistler. |
|
图 18 AUC盒形图 Fig. 18 AUC box |
(1) Gray_Blue_L特征分类器在精度(Precision)、召回率(Recall)和F1score三个指标上的表现最好且识别性能最稳定.比如图 17a最左边的精度图，其横轴是不同的图像特征提取方法，纵轴是精度，仔细观察发现：Gray_Blur_L特征对应的箱体的高度最小，说明该特征分类器的精度最稳定；同时，Gray_Blur_L箱体的中位线均高于Gray_Blur特征和Gray特征的箱体的中位线，表明Gray_Blur_L特征分类器在识别闪电哨声波时，其精度最高；按照同样的方法观察图 17a中间的图和最右边的图可得结论：在识别闪电哨声波时，Gray_Blur_L特征分类器的召回率和F1score指标的值最优且分类性能最稳定；按照同样的方式观察如图 17b的三幅子图，可以得到如下结论：在对非闪电哨声波进行分类时，Gray_Blue_L特征分类器在精度、召回率和F1score值三个指标的值最高，且性能最稳定.
(2) 观察三个分类器的10000个AUC指标的盒形图，如图 18所示：其横轴是不同的特征提取方法，纵轴是AUC值.发现：Gray_Blue_L的箱体高度最小中位线的位置最高.这说明Gray_Blue_L特征所训练的分类器不仅仅波动较小且具有最好的排序能力.
总之，由Gray_Blur_L图像特征学习得到的分类器，不仅仅具有更可靠的分类效果且分类性能更稳定.
3.4.4 识别效果显著性差异评价为了检验不同的特征分类器的性能是否存在明显差异，采用T检验方法对显著性差异进行定量评价，其显著性P值越高，表明显著差异性越小，通常阈值是0.05，其含义是若差异性小于0.05则认为存在明显差异；若大于0.05则认为两组实验不存在明显差异.其结果如表 5—8所示.
|
|
表 5 T检验P值(Precision) Table 5 T tests for Precision |
|
|
表 6 T检验P值(Recall) Table 6 T tests for Recall |
|
|
表 7 T检验P值(F1Score) Table 7 T tests for F1Score |
|
|
表 8 T检验P值(AUC) Table 8 T tests for AUC |
首先，对三种不同特征分类器的精度(Precision)指标的分布数据，两两进行T检验，得到如表 5的检验值.比如对Gray和Gray_Blur的数据进行T检验，得到其P值是0，表明：Gray特征分类器和Gray_Blur特征分类器在分类精度方面存在明显差异.再比如，对Gray特征分类器的自身进行比较检验，即对Gray特征分类器和Gray特征分类器进行T检验的P值是1，表明：数据之间无显著性差异.然后，按照上述方法观察表 6—8，发现：三种不同特征的分类器在召回率方面、排序能力等方面均存在明显差异.
总之，采用Gray特征，Gray_Blur特征，Gray_Blur_L特征进行闪电哨声波识别，其识别效果在精度、召回率、F1值和AUC值均存在明显差异，且Gray_Blur_L特征分类器的分类效果最佳.
4 讨论上述实验表明基于张衡卫星SCM数据开展自动化识别闪电哨声波的算法研究具有一定的效果.算法方案中的模糊卷积核和L特征卷积核对闪电哨声波识别具有非常重要的影响.本章节将对其产生的影响开展较深入的讨论和分析，同时对闪电哨声波自动识别的延伸研究进行展望.
4.1 模糊卷积核实验结果表明采用模糊卷积核提取图像特征进行闪电分类的效果优于采用原始频谱图为特征的识别效果，本小节将从视觉观察角度和特征可分角度对其展开讨论.
4.1.1 视觉观察角度原始频谱图具有明显的垂直信息痕迹，如图 19所示：大量的垂直信息的存在对L形态特征的识别产生了巨大的影响，而图 20是采用模糊卷积核对垂直信息进行模糊处理后的图像，其垂直信息被减弱，同时凸显了其L形态特征.
|
图 19 含有闪电信息的原始灰度图 Fig. 19 Original gray image containing lightning whistlers |
|
图 20 采用模糊卷积核处理后的图像 Fig. 20 The resulting image processed by fuzzy convolution kernel |
为了探索模糊卷积处理后的图像特征与原始灰度图的特征是否具有不同的可分性.本小节将采用主成分分析(Principal Components Analysis, PCA)的降维技术将316个闪电哨声波样本和1000个随机的非闪电哨声波样本降至2维，并绘制出其在2维空间的分布，以观察特征的可分性趋势，如图 21所示.
|
图 21 特征分布趋势图 Fig. 21 Character distribution |
图 21a是闪电哨声波与非闪电哨声波数据的原始特征分布，其横坐标和纵坐标表示降维后的2个维度；图 21b是被模糊卷积处理后的特征分布；其中1代表闪电哨声波，-1是非闪电哨声波，观察图 21发现：原始图像特征有大量的正负样本聚集在一起的情况，尤其是图 21a中的圆形区域，这些叠加在一起的正负样本增加了识别闪电哨声波和非闪电哨声波的难度，导致分类器识别闪电的性能下降；相比之下，图 21b的圆形区域，正负样本聚集在一起的数量较前者有大幅减少，该现象使得正负样本的区分提高，增加了闪电哨声波和非闪电哨声波图像特征的区分度，使得分类器识别闪电哨声波的能力增强.
4.2 L特征卷积核对闪电哨声波识别的影响从实验结果发现在模糊卷积处理图像的基础上，再采用9×9的L特征卷积核处理，使得闪电哨声波分类的效果更优，接下来针对L特征卷积核的识别效果和尺度选择从两个方面进行深入讨论.
4.2.1 视觉观察角度从图 22中发现相比较原始图像，模糊核提取的L特征更为明显，如图 22a和图 22b所示；进一步观察，采用3×3的L卷积核提取的L特征如图 22c所示，由于太暗，将图 22c的矩形区域进行放大如图 22j所示，发现大部分的细节信息被滤除，隐隐约约可看到L特征(如图 22j的椭圆区域所示)，说明其只保留了L特征的主轮廓；随着L卷积核尺度变大，L特征周围的细节信息越来越多，如图 22(d，e，f)所示；随着L卷积核的尺度进一步增大，不仅L特征周围的细节强度越来越明显，且与L区域无关的地方亦出现了较多的细节.
|
图 22 不同卷积核提取的图像特征(含有闪电哨声波的图像) (a) 原始灰度图; (b) 模糊卷积核处理后的图像; (c) 尺度3×3的L卷积核处理b后的图像；(d) 尺度5×5的L卷积核处理b后的图像；(e) 尺度7×7的L卷积核处理b后的图像; (f) 尺度9×9的L卷积核处理b后的图像; (g) 尺度11×11的L卷积核处理b后的图像; (h) 尺度13×13的L卷积核处理b后的图像; (i) 尺度15×15的L卷积核处理b后的图像. Fig. 22 Image features extracted by different convolution kernel (containing lightning whistlers) (a) The original grayscale image; (b) The image processed by the fuzzy convolution kernel; (c) The image obtained by processing the (b) with 3×3 L convolution kernel; (d) The image obtained by processing the (b) with 5×5 L convolution kernel; (e) The image obtained by processing the (b) with 7×7 L convolution kernel; (f) The image obtained by processing the (b) with 9×9 L convolution kernel; (g) The image obtained by processing the (b) with 11×11 L convolution kernel; (h) The image obtained by processing the (b) with 13×13 L convolution kernel; (i) The image obtained by processing the (b) with 15×15 L convolution kernel. |
进一步观察图 23，小尺度的L卷积核能够滤除细节信息，如图 23(c，d，e)所示；随着L卷积核尺度的不断增加，其细节信息越来越明显，如图 23(f，g)所示；当尺度进一步增加时，放大的干扰信息产生了与L特征相似的形态，如图 23i的圆形区域，即大尺度L卷积核容易将该区域的噪声放大，容易影响识别性能.因此，小尺度的L卷积核滤除L特征周围太多的细节信息，太大的L卷积核又容易放大非闪电图像中的噪声，如何自动选择尺度合适的L卷积核是一个需要进一步思考的问题.
|
图 23 不同卷积核提取的图像特征(含有非闪电哨声波的图像) (a) 原始灰度图；(b) 模糊卷积核处理后的图像; (c) 尺度3×3的L卷积核处理b后的图像；(d) 尺度5×5的L卷积核处理b后的图像；(e) 尺度7×7的L卷积核处理b后的图像; (f) 尺度9×9的L卷积核处理b后的图像; (g) 尺度11×11的L卷积核处理b后的图像; (h) 尺度13×13的L卷积核处理b后的图像; (i) 尺度15×15的L卷积核处理b后的图像. Fig. 23 Image features extracted by different convolution kernel (not containing lightning whistlers) (a) The original grayscale image; (b) The image processed by the fuzzy convolution kernel; (c) The image obtained by processing the (b) with 3×3 L convolution kernel; (d) The image obtained by processing the (b) with 5×5 L convolution kernel; (e) The image obtained by processing the (b) with 7×7 L convolution kernel; (f) The image obtained by processing the (b) with 9×9 L convolution kernel; (g) The image obtained by processing the (b) with 11×11 L convolution kernel; (h) The image obtained by processing the (b) with 13×13 L convolution kernel; (i) The image obtained by processing the (b) with 15×15 L convolution kernel. |
为了定性地比较不同特征的可分性，首先采用不同的特征提取方式提取316个闪电哨声波样本和1000个随机的非闪电哨声波样本的图像特征，并采用PCA技术将特征维度降至二维，最后，将其绘制在二维空间，以观察其相应的特征的分布情况，并分析其可分性的趋势，其中+1表示闪电样本数据，-1表示非闪电样本数据，如图 24所示.
|
图 24 图像特征降至二维的示意图 (a) 模糊卷积核处理后的图像; (b)尺度3×3的L卷积核处理a后的图像；(c) 尺度5×5的L卷积核处理a后的图像；(d) 尺度7×7的L卷积核处理a后的图像; (e) 尺度9×9的L卷积核处理a后的图像; (f) 尺度11×11的L卷积核处理a后的图像; (g) 尺度13×13的L卷积核处理a后的图像; (h) 尺度15×15的L卷积核处理a后的图像. Fig. 24 Represent the image characters at two-dimensional space (a) The images obtained by processing the original image with fuzzy convolution kernel; (b) The images obtained by processing (a) with 3×3 L convolution kernel; (c) The images obtained by processing (a) with 5×5 L convolution kernel; (d) The images obtained by processing (a) with 7×7 L convolution kernel; (e) The images obtained by processing (a) with 9×9 L convolution kernel; (f) The images obtained by processing (a) with 11×11L convolution kernel; (g) The images obtained by processing (a) with 13×13 L convolution kernel; (h) The images obtained by processing (a) with 15×15 L convolution kernel. |
通过图 24可发现：模糊卷积核提取的图像特征如图 24a所示，虽然整体上看大部分的正负样本均已分开，但在中间仍然存在正负样本叠加在一起的情况；图 24b是在通过先采用模糊卷积处理图像，再采用3×3的L卷积核处理图像得到的特征在二维空间的分布效果，与图 24a相比，正负样本之间的距离拉开了，随着L卷积核尺度的增加，正负类之间的可分性均优于图 24a，但与图 24b相比其差异不大，如图 24(c，d，f)所示；然而值得一提的是，闪电哨声波样本的类内越来越积聚，表明闪电哨声波样本的类内差越来越小，该情况有利于闪电事件的识别，如图 24f；随着L特征卷积核的尺度的进一步增大，发现部分闪电哨声波样本和非闪电样哨声波样本存在相互重叠的现象，如图 24(g，h)中的圆形区域所示.而且随着尺度的增加，这个重叠现象越来越严重，如图 24h所示.此外，其非闪电哨声波数据的类内越来越疏散，表明非闪电哨声波样本的类内差越来越大，该情况不利于对非闪电哨声波进行分类.
总之，尺度合适的L卷积核能够较好地对闪电哨声波与非闪电哨声波进行分类.太小的尺度虽然能够提高分类效果，但是闪电哨声波的类内差大，不利于闪电哨声波的识别；太大的尺度容易造成非闪电哨声波和闪电哨声波在特征空间发生重叠，导致类间方差越来越小；亦导致非闪电哨声波类内方差变大；类间方差小和类内方差大将严重影响分类器的分类效果.
4.2.3 统计盒形图的角度利用不同尺度的L卷积核提取的图像特征并训练闪电哨声波分类器，计算反映分类器排序能力的AUC值，重复实验10000次，得到该特征下的分类器的AUC指标的分布情况，如图 25所示：其横轴是不同的特征提取方法，最左边是采用模糊卷积核提取的图像特征，接下来依次是3×3 L形态卷积核、5×5 L形态卷积核、…、15×15 L形态卷积核.发现原始灰度特征的AUC的箱体远远高于基于L形态卷积核的箱体.尤其在卷积核是13×13的时候，其AUC箱体不仅仅最小而且其中位数最高.说明采用尺度为13×13的L卷积核提取的特征训练得到的分类器，具有最稳定最优的排序能力.
|
图 25 不同特征训练的分类器的AUC Fig. 25 Use different features to train the different classifier |
同样，采用不同尺度的L卷积核提取的特征训练分类器，重复10000次，可分别获得10000个精度指标值(Precision)、召回率值和F1score值，其数据分布如图 26a所示，图中Blur表示采用模糊卷积核提取的特征，L3表示3×3的L形卷积核，以此类推，L15表示15×15的L形卷积核.观察发现：在识别闪电哨声波的过程中，未采用L卷积核提取特征的分类器，其Recall指标的箱体均高于采用L形卷积特征的箱体，特别是在图 26a第2幅图的L9的位置，该位置意味着卷积核的尺度是9×9，同时发现该位置的中位线均高于其他位置的箱体的中位线，表明采用该尺度卷积核提取的特征训练分类器，其性能在召回率上不仅仅波动最小而且精度最优.同时发现，Recall数据的箱体均在L3和L15位置，箱体高度变高，且中位线下降.意味着尺度太小的L形卷积核，如尺度为3×3，和尺度太大的卷积核，如尺度为15×15的L形卷积核，均使得闪电哨声波识别的性能下降.上述同样的性能亦体现在F1score的分布图上.但其精度变化不是很大.图 26b呈现的是识别非闪电哨声波的三个指标的数据，情况分别是：在精度指标上，L9位置的中位线达到最高，表明该尺度的卷积核效果在该位置上最佳，与该位置相比，最小的尺度(L3)和最大的尺度(L15)效果均不理想；在召回率指标上，卷积核尺度对其影响不大；在F1score，最小尺度的卷积核效果较差.
|
图 26 不同特征训练的分类器的三个指标分布情况 (a) 识别闪电哨声波；(b) 识别非闪电哨声波. Fig. 26 Performance of the classifiers (a) Lightning whistlers; (b) Without lightning whistlers. |
总之，基于卷积核提取的图像特征不仅仅具有更好的分类性能且比原始灰度的波动小，分类性更稳定.
4.3 基于电磁卫星SCM数据的闪电哨声波识别研究的延伸讨论基于电磁卫星的SCM数据的闪电哨声波识别方案可行且鲁棒性较好，以此可以获得大量的闪电哨声波数据，在此基础上，可开展下列四个方面的研究：
(1) 通过式(13)得到质子回旋频率(fc)：
|
(13) |
其中q是质子电荷，为1.60217733×10-19C，m是质子的质量，为1.672621637×10-27kg; B是磁感应强度，单位是T.绘制质子回旋频率，如图 27中箭头所指.以质子回旋频率作为参考，分析闪电哨声波相对于质子回旋频率的位置，借此，进一步分析全球闪电哨声波随地点、季节、年等时空变化的规律.
|
图 27 质子回旋频率 Fig. 27 Proton cyclotron frequency |
(2) 从图 27发现闪电哨声波存在许多不同类型的色散形态.该现象与闪电哨声传播环境有关.闪电发射的脉冲辐射沿着地球-电离层波导传播.在电离层底层一部分辐射溢出地球-电离层波导传到电离层，并以哨声模传播.因为电子数陡峭的正梯度，向上传播的哨声模波动方向角在穿越电离层时被折射，在F层顶端变成近似垂直.在穿越电离层时，波动功率谱密度因为库伦损失被减弱.当波动到达等离子体层的时候，最初的VLF辐射已经远离它的源，由等离子体层电子密度梯度和磁场梯度引导传播(曹晋滨等, 2009).当上述环境发生变化时，其呈现的色散形态亦不同.可通过对闪电哨声波色散形态进一步开展聚类分析来挖掘全球电离层环境的特点.
(3) Fiser等(2010)提出哨声波幅度可视作经过某固定距离闪电的电流强度的函数，并指出欧洲的闪电哨声波的幅度与闪电源和卫星轨道的位置有很大关系.可以结合地面闪电观测系统与卫星观测电场数据，探索闪电源到卫星之间的距离与闪电哨声波幅度的关系.
(4) 本方案可进一步拓展到从卫星数据中去识别其他具有明确形态特征的空间物理现象.比如嘶声、合声、地面甚低频发射机VLF波动、地面电力系统的谐频辐射、电磁离子回旋波等.这些空间现象的频谱具有明显的形态特征，自动识别这些非震的空间现象有助于剥离非震的电磁扰动，推动自动识别地震电磁扰动的发展.
(5) 本文介绍的方法是提取如图 28所示的时频图像中的L形态特征(如图 29所示)并进一步识别该时频图中是否存在闪电哨声波，未直接对闪电哨声波进行描迹.基于闪电哨声波轨迹提取的截止频率、持续时间、色散阶数等信息对认知哨声波的传播过程、获取背景参数信息至关重要.未来将以本文提出的方法为基础，实现基于时频图的闪电哨声波描迹：当识别某时频图中存在闪电哨声波后，首先将如图 28所示的彩色时频图像中的每个像素点从RGB颜色空间转化到HSV颜色空间，设置色调阈值T0；再通过本文提出的方法将彩色时频图像转换成L特征图，如图 29所示，并设置阈值T1；若某像素点在HSV颜色空间中的H值高于该阈值T0，同时该点在L特征图中像素值高于T1，则表明该点为闪电哨声波轨迹点，以完成自动描迹功能.
|
图 28 含有闪电哨声波的时频图 Fig. 28 Time-Frequency profile |
本文在张衡一号卫星的SCM数据上开展自动化识别闪电哨声波的研究.根据频谱数据的特点和闪电哨声波的L形态特征设计了模糊卷积核和L形态卷积核，进一步增强了频谱图中的闪电形态特征，并采用SVM分类器对特征进行分类，使得其识别性能在精度、召回率以及排序能力等方面均处于94%以上.
然而，由于闪电哨声波样本(正样本)量太少，非闪电哨声波样本(负样本)量大，为了避免分类器在分类的过程中容易出现偏好负样本的情况，文中采用的是降采样负样本的策略.为了得到更为鲁棒的闪电哨声波识别效果，需要增加正样本数量，以平衡正负样本量.因此，未来的工作主要集中在三个方面：
(1) 在本次探索基础上进一步增加正样本量，并探索采用卷积神经网络(CNN)的自动化特征提取方法；
(2) 当获得大量的闪电哨声波事件后进一步提取其频率和幅度的变化并联合其他载荷数据，建立闪电哨声波频率/幅度与等离子体离子密度、离子温度、离子速度、等离子体电子温度、电子密度之间的数学模型；
(3) 本研究提出的识别方案，可进一步拓展应用于从张衡卫星的数据中自动识别人工发射源、磁暴等空间事件；
(4) 鉴于我国地面ELF/VLF接收机系统获得的突破性进展(Chen et al., 2016, 2017; Zhou et al., 2020; Yi et al., 2020)，下一步将利用张衡一号卫星数据和我国地面的ELF/VLF接收机观测数据开展闪电哨声波联合观测研究.
致谢 本文特别感谢来自应急管理部国家自然灾害防治研究院的张衡一号卫星团队的所有成员为本研究提供的技术服务支持.
Ali Ahmad U, Kasahara Y, Matsuda S, et al. 2019. Automatic detection of lightning whistlers observed by the plasma wave experiment onboard the Arase satellite using the OpenCV library. Remote Sensing, 11(15): 1785. DOI:10.3390/rs11151785 |
Bayupati I P A, Kasahara Y, Goto Y. 2012. Study of dispersion of lightning whistlers observed by Akebono satellite in the earth's plasmasphere. IEICE Transactions on Communications, E95.B(11): 3472-3479. DOI:10.1587/transcom.E95.B.3472 |
Bresenham J. 1977. A linear algorithm for incremental digital display of circular arcs. Communications of the ACM, 20(2): 100-106. DOI:10.1145/359423.359432 |
Christian, Hugh J. 2003. Global frequency and distribution of lightning as observed from space by the Optical Transient Detector. Journal of Geophysical Research Atmospheres, 108(D1): ACL 4-1-ACL 4-15. |
Cai J T, Zhao G Z, Zhan Y, et al. 2007. The study on ionospheric disturbances during earthquakes. Progress in Geophysics (in Chinese), 22(3): 695-701. DOI:10.3969/j.issn.1004-2903.2007.03.006 |
Cao J B, Yan C X, Lu L, et al. 2009. Non-seismic induced electromagnetic waves in the near earth space. Earthquake (in Chinese), 29(S1): 17-25. |
Chen Y P, Yang G B, Ni B B, et al. 2016. Development of ground-based ELF/VLF receiver system in Wuhan and its first results. Advances in Space Research, 57(9): 1871-1880. DOI:10.1016/j.asr.2016.01.023 |
Chen Y P, Ni B B, Gu X D, et al. 2017. First observations of low latitude whistlers using WHU ELF/VLF receiver system. Science China Technological Sciences, 60(1): 166-174. DOI:10.1007/s11431-016-6103-5 |
Chmyrev V M, Isaev N V, Bilichenko S V, et al. 1989. Observation by space-borne detectors of electric fields and hydromagnetic waves in the ionosphere over an earthquake centre. Physics of the Earth and Planetary Interiors, 57(1-2): 110-114. DOI:10.1016/0031-9201(89)90220-3 |
Clilverd M A, Nunn D, Lev-Tov S J, et al. 2002. Determining the size of lightning-induced electron precipitation patches. Journal of Geophysical Research: Space Physics, 107(A8): SIA 10-1-SIA 10-11. DOI:10.1029/2001JA000301 |
Fiser J, Chum J, Diendorfer G, et al. 2010. Whistler intensities above thunderstorms. Annales Geophysicae, 28(1): 37-46. DOI:10.5194/angeo-28-37-2010 |
Gurnett D A, Kurth W S, Cairns I H, et al. 1990. Whistlers in Neptune's magnetosphere: Evidence of atmospheric lightning. Journal of Geophysical Research: Space Physics, 95(A12): 20967-20976. DOI:10.1029/JA095iA12p20967 |
Helliwell R A. 1965. Whistlers and Related Ionospheric Phenomena. Stanford, California: Stanford University Press.
|
Larkina V I, Nalivaiko A V, Gershenzon N I, et al. 1983. Intercosmos-19 observations of VLF emissions associated with seismic activity. Geomagnetizm i Aeronomiia, 23: 842-846. |
Li L Y, Yang J Y, Cao J B, et al. 2011. Statistical backgrounds of topside-ionospheric electron density and temperature and their variations during geomagnetic activity. Chinese Journal of Geophysics (in Chinese), 54(10): 2437-2444. DOI:10.3969/j.issn.0001-5733.2011.10.001 |
Li L Y, Cao J B, Yang J Y, et al. 2015. Semiannual and solar activity variations of daytime plasma observed by DEMETER in the ionosphere-plasmasphere transition region. Journal of Geophysical Research: Space Physics, 120(12): 10640-10653. DOI:10.1002/2015JA021102 |
Liao L, Zhao S F, Shen X H, et al. 2019. Characteristic analysis and full wave simulation of electrical field for China seismo-electromagnetic satellite observations radiated from VLF transmitter. Chinese Journal of Geophysics (in Chinese), 62(4): 1210-1217. DOI:10.6038/cjg2019M0504 |
Liu Y M, Wang J S, Xiao Z, et al. 2006. A possible mechanism of typhoon effects on the Ionospheric F2 layer. Chinese Space Science (in Chinese), 26(2): 92-97. |
Molchanov O A, Mazhaeva O A, Golyavin A N, et al. 1993. Observation by the Intercosmos-24 satellite of ELF-VLF electromagnetic emissions associated with earthquakes. Annales Geophysicae, 11(5): 431-440. |
Nilupar M, Zhang Y X. 2012. Recent research advances on ionospheric disturbances from seismo-electromagnetic satellite. Earthquake (in Chinese), 32(1): 103-117. |
Oike Y, Kasahara Y, Goto Y. 2014. Spatial distribution and temporal variations of occurrence frequency of lightning whistlers observed by VLF/WBA onboard Akebono. Radio Science, 49(9): 753-764. DOI:10.1002/2014RS005523 |
Parrot M, Mogilevsky M M. 1989. VLF emissions associated with earthquakes and observed in the ionosphere and the magnetosphere. Physics of the Earth and Planetary Interiors, 57(1-2): 86-99. DOI:10.1016/0031-9201(89)90218-5 |
Parrot M, Berthelier J J, Blecki J, et al. 2015. Unexpected Very Low Frequency (VLF) Radio Events Recorded by the Ionospheric Satellite DEMETER. Surveys in Geophysics, 36(3): 483-511. DOI:10.1007/s10712-015-9315-5 |
Parrot M, Pinçon J L, Shklyar D. 2019. Short-fractional hop whistler rate observed by the low-altitude satellite DEMETER at the end of the solar cycle 23. Journal of Geophysical Research: Space Physics, 124(5): 3522-3531. DOI:10.1029/2018JA026176 |
Pulinets S A, Legen'ka A D. 2003. Spatial-temporal characteristics of large scale disturbances of electron density observed in the Ionospheric F-region before strong earthquakes. Cosmic Research, 41(3): 221-230. DOI:10.1023/A:1024046814173 |
Sánchez A D V. 2003. Advanced support vector machines and kernel methods. Neurocomputing, 55(1-2): 5-20. DOI:10.1016/S0925-2312(03)00373-4 |
Santolík O, Parrot M, Inan U S, et al. 2009. Propagation of unducted whistlers from their source lightning: A case study. Journal of Geophysical Research: Space Physics, 114(A3): A03212. DOI:10.1029/2008JA013776 |
Shen X H, Zong Q-G, Zhang X M. 2018. Introduction to special section on the China Seismo-Electromagnetic Satellite and initial results. Earth Planet. Phys., 2(6): 439-443. DOI:10.26464/epp2018041 |
Tian M, Xu J S, Ma S Y. 1991. KDX-87 Receiving and analysing system for VLF wide-band direction finding. Acta Geophysica Sinica (in Chinese), 34(1): 11-19. |
Wang Q, Huang J P, Zhang X M, et al. 2018. China Seismo-Electromagnetic Satellite search coil magnetometer data and initial results. Earth Planet. Phys., 2(6): 462-468. DOI:10.26464/epp2018044 |
Wang T Z, Liu Z H, Yan C G. 1982. A preliminary analysis on whistler observation over Zhanjiang region, China. Acta Geophysica Sinica (in Chinese), 25(5): 391-397. |
Xia L D, Wang Y S. 1994. A statistical analysis of the observations of whistler at Sanya Area. Journal of China University of Science and Technology (in Chinese), 24(3): 368-372. |
Xiao S G, Zhang D H, Xiao Z. 2007. Study on the Detectability of Typhoon-Generated Acoustic-Gravity Waves. Chinese Journal of Space Science (in Chinese), 27(1): 35-40. |
Xu J S, Tian M, Ma S Y, et al. 1989. Multi-stationed wideband direction finding measurements for whistlers at geomagnetic latitudes below 20° in China and some early results. Acta Geophysica Sinica (in Chinese), 32(2): 125-134. |
Yan R, Shen X H, Huang J P, et al. 2018. Examples of unusual ionospheric observations by the CSES prior to earthquakes. Earth Planet. Phys.,, 2(6): 515-526. DOI:10.26464/epp2018050 |
Yang Y Y, Zhima Z R, Shen X H, et al. 2020. The first intense geomagnetic storm event recorded by the China Seismo-Electromagnetic Satellite. Space Weather, 18(1): e2019SW002243. DOI:10.1029/2019SW002243 |
Yi J, Gu X D, Cheng W, et al. 2020. A detailed investigation of low latitude tweek atmospherics observed by the WHU ELF/VLF receiver: 2. Occurrence features and associated ionospheric parameters. Earth Planet. Phys., 4(3): 238-245. DOI:10.26464/epp2020023 |
Yuan S G, Zhu X H, Huang J P. 2018. System design and key technology of China Seismo-Electromagnetic Satellite. Journal of Remote Sensing (in Chinese), 22(S1): 32-38. |
Záhlava J, Němec F, Santolík O, et al. 2018. Longitudinal dependence of whistler mode electromagnetic waves in the earth's inner magnetosphere. Journal of Geophysical Research: Space Physics, 123(8): 6562-6575. DOI:10.1029/2018JA025284 |
Zhang Y J. 2016. A Course of Image Processing and Analysis (in Chinese). 2nd ed. Beijing: People's Posts and Telecommunications Press: 242.
|
Zhima Z R, Huang J P, Shen X H, et al. 2020. Simultaneous Observations of ELF/VLF Rising-Tone Quasiperiodic Waves and Energetic Electron Precipitations in the High-Latitude Upper Ionosphere. Journal of Geophysical Research: Space Physics, 125(5): 1-21. DOI:10.1029/2019ja027574 |
Zhao S F, Zhou C, Shen X H, et al. 2019. Investigation of VLF transmitter signals in the ionosphere by ZH-1 observations and full-wave simulation. Journal of Geophysical Research: Space Physics, 124(6): 4697-4709. DOI:10.1029/2019JA026593 |
Zhou B, Yang Y Y, Zhang Y T, et al. 2018. Magnetic field dataprocessing methods of the China Seismo-Electromagnetic Satellite. Earth Planet. Phys., 2(6): 455-461. DOI:10.26464/epp2018043 |
Zhou R X, Gu X D, Yang K X, et al. 2020. A detailed investigation of low latitude tweek atmospherics observed by the WHU ELF/VLF receiver: I. Automatic detection and analysis method. Earth Planet. Phys., 4(2): 120-130. DOI:10.26464/epp2020018 |
蔡军涛, 赵国泽, 詹艳, 等. 2007. 地震期间电离层扰动现象研究. 地球物理学进展, 22(3): 695-701. DOI:10.3969/j.issn.1004-2903.2007.03.006 |
曹晋滨, 燕春晓, 路立, 等. 2009. 地球近地空间非震电磁扰动. 地震, 29(S1): 17-25. |
李柳元, 杨俊英, 曹晋滨, 等. 2011. 顶部电离层电子密度和温度的统计背景及其地磁活动变化. 地球物理学报, 54(10): 2437-2444. DOI:10.3969/j.issn.0001-5733.2011.10.001 |
廖力, 赵庶凡, 申旭辉, 等. 2019. ZH-1卫星观测的VLF人工源信号特征分析与全波模拟. 地球物理学报, 62(4): 1210-1217. DOI:10.6038/cjg2019M0504 |
刘依谋, 王劲松, 肖佐, 等. 2006. 台风影响电离层F2区的一种可能机制. 空间科学学报, 26(2): 92-97. DOI:10.3969/j.issn.0254-6124.2006.02.003 |
尼鲁帕尔·买买吐孙, 张永仙. 2012. 地震电磁卫星电离层扰动研究进展综述. 地震, 32(1): 103-117. DOI:10.3969/j.issn.1000-3274.2012.01.010 |
田茂, 徐继生, 马淑英. 1991. KDX-87型甚低频宽带定向接收与微机分析系统. 地球物理学报, 34(1): 11-19. DOI:10.3321/j.issn:0001-5733.1991.01.002 |
王挺柱, 刘曾厚, 严成高. 1982. 湛江地区哨声观测的初步分析. 地球物理学报, 25(5): 391-397. DOI:10.3321/j.issn:0001-5733.1982.05.001 |
夏利东, 王友善. 1994. 三亚地区历年哨声观测资料的统计分析. 中国科学技术大学学报, 24(3): 368-372. |
肖赛冠, 张东和, 肖佐. 2007. 台风激发的声重力波的可探测性研究. 空间科学学报, 27(1): 35-40. DOI:10.3969/j.issn.0254-6124.2007.01.007 |
徐继生, 田茂, 马淑英, 等. 1989. 磁纬20°以下地区哨声多台宽带定向观测及其初步结果. 地球物理学报, 32(2): 125-134. |
袁仕耿, 朱兴鸿, 黄建平. 2018. 电磁监测试验卫星(张衡一号)系统设计与关键技术. 遥感学报, 22(S1): 32-38. |
章毓晋. 2016. 图像处理和分析教程. 2版. 北京: 人民邮电出版社: 242.
|