张衡一号感应磁力仪数据闪电哨声波自动识别

袁静， 王桥， 杨德贺， 刘芹芹， 泽仁志玛， 申旭辉

1 防灾科技学院， 河北三河 065421 2 应急管理部国家自然灾害防治研究院， 北京 100085

0 引言

利用电磁卫星监测地震始于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的数据量.如何从如此庞大的观测数据中自动识别闪电哨声波事件尤为关键和紧迫.

图1 AKEBONO卫星闪电形态图例(Bayupati et al., 2012)Fig.1 Lightning whistlers of AKEBONO (Bayupati et al., 2012)

目前基于海量闪电哨声波数据开展全球时空分布规律和相关参数的研究较少，主要存在以下两个方面的原因：一是需要足够的闪电哨声波事件；二是缺少有效的闪电哨声波自动识别算法，依靠人工从大量的电磁观测数据中获取闪电哨声波事件是一件极具有挑战性的工作.比如：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所示.

图2 DEMETER卫星闪电形态图例 (Parrot et al., 2015, 2019)Fig.2 Lightning whistlers of DEMETER (Parrot et al., 2015, 2019)

表1 闪电哨声波分类(Helliwell, 1965)Table 1 Type and principal definitions of lightning whistlers (Helliwell, 1965)

目前基于电磁卫星闪电哨声波自动识别的研究主要是针对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.

2 闪电哨声波自动识别算法

闪电哨声波识别算法的目标是识别每个时频图中是否有闪电哨声波.目前深度学习模型是机器学习的主流模型，但仍需要大量的数据样本作为基础.而本文目前收集的闪电哨声波数据量较少，基于此仍然用传统的识别算法技术探索闪电哨声波识别方案，方案包括训练过程和识别过程，如图8所示：实线是训练流程，虚线是识别流程.训练流程的主要目的是获得闪电哨声波的识别模型，识别流程的主要目的是应用闪电哨声波识别模型.训练过程包括：灰度化处理、尺度处理，模糊卷积处理，L形态卷积处理，训练SVM模型.其识别流程包括：灰度化处理、尺度处理，模糊卷积处理，L形态卷积处理，应用SVM模型进行闪电识别.

图8 识别算法(实线是模型训练流；虚线是模型识别流)Fig.8 Recognition algorithm (the solid line represents training process; the dotted line denotes recognition process)

2.1 灰度化处理

观察原始的样本频谱图，如图9所示，发现识别闪电哨声波的依据是频谱图中是否存在L形态特征，无需颜色特征，则首先采用灰度化处理消除颜色的影响.对图9按照公式(3)进行灰度化处理(章毓晋, 2016)，得到如图10所示的灰度图.

图9 彩色频谱图(有闪电)Fig.9 Colorful spectrum (lightning whistlers)

图10 灰度化(有闪电)Fig.10 Grayscale (lightning whistlers)

Gray=RGB.R×0.3+RGB.G×0.59

+RGB.B×0.11，

(3)

其中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

2.3 模糊卷积

观察图12发现该图像条纹状纹理比较严重，表明该图具有明显的阶跃边缘.这些边缘信息不利于L形态特征的识别.为了弱化阶跃边缘对识别的影响，本文设计如图13a所示的模糊卷积核，用其对图12进行卷积处理，得到如图13b所示的效果：其阶跃边缘信息被弱化了.

图13 模糊卷积(a) 模糊卷积核； (b) 模糊卷积效果.Fig.13 Fuzzy convolution(a) Kernel； (b) Result.

2.4 L特征提取

为了进一步增强闪电哨声波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.

2.5 SVM分类

将图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)所示：

wx+b=0，

(4)

其中x是在超平面上的样本点，b是偏置.

若存在使得分类间隔最大的超平面，其满足下列条件：闪电哨声波的样本正确地落在闪电哨声波的分类区，则有公式(5)；

(wxi+b)≥+1，yi=1.

(5)

非闪电哨声波的样本正确落在非闪电哨声波的分类区，则有公式(6)；

(wxi+b)≤-1，yi=-1，

(6)

公式(5)和公式(6)合成一个公式(7)：

yi(wxi+b)-1≥0.

(7)

SVM将在众多的分类超平面中，自动寻找使得特征空间间隔最大的超平面.其数学模型表示为公式(8)：

max 1/‖w‖ s.t.yi(wxi+b)-1≥0

(i=1，2，…,n).

(8)

训练模型的目的是找到满足式(8)的w.

由于大部分情况下，数据无法满足线性可分的要求，此时满足上述条件的超平面不存在，对于这种非线性可分的问题，SVM 提供核函数k(xi，xj)代替公式(5)，将原数据映射到高维空间，使得其在高维空间线性可分(Sánchez and David, 2003).本文中的核函数采用的是多项式核函数如式(9)所示：

(9)

3 实验和分析

3.1 实验方案

由于闪电哨声波的样本与非闪电哨声波的样本不均衡，按照传统的机器学习的惯例：正负样本的比例不能超过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)所示：

Precision =TP/(FP+TP)，

(10)

含义：在预测出的所有正例中(TP+FP)，真正例(TP)占的比例，通常该指标越大越好.

(2)召回率(Recall)的定义如式(11)所示：

Recall =TP/(TP+FN) ，

(11)

含义：在数据集中所有的正例中(TP+FN)，真正例(TP)占的比例，通常该指标越大越好.

(3)F1值(F1-Score)的定义如式(12)所示：

F1=2/(1/ Precision +1/Recall) ，

(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.

3.4.2 总体识别精度评价

本小节主要针对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

通过观察表3和表4可发现：

(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所示.

首先，对三种不同特征分类器的精度(Precision)指标的分布数据，两两进行T检验，得到如表5的检验值.比如对Gray和Gray_Blur的数据进行T检验，得到其P值是0，表明：Gray特征分类器和Gray_Blur特征分类器在分类精度方面存在明显差异.再比如， 对Gray特征分类器的自身进行比较检验，即对Gray特征分类器和Gray特征分类器进行T检验的P值是1，表明：数据之间无显著性差异.然后，按照上述方法观察表6—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

总之，采用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

4.1.2 特征可分的角度

为了探索模糊卷积处理后的图像特征与原始灰度图的特征是否具有不同的可分性.本小节将采用主成分分析(Principal Components Analysis, PCA)的降维技术将316个闪电哨声波样本和1000个随机的非闪电哨声波样本降至2维，并绘制出其在2维空间的分布，以观察特征的可分性趋势，如图21所示.

图21a是闪电哨声波与非闪电哨声波数据的原始特征分布，其横坐标和纵坐标表示降维后的2个维度；图21b是被模糊卷积处理后的特征分布；其中1代表闪电哨声波，-1是非闪电哨声波，观察图21发现：原始图像特征有大量的正负样本聚集在一起的情况，尤其是图21a中的圆形区域，这些叠加在一起的正负样本增加了识别闪电哨声波和非闪电哨声波的难度，导致分类器识别闪电的性能下降；相比之下，图21b的圆形区域，正负样本聚集在一起的数量较前者有大幅减少，该现象使得正负样本的区分提高，增加了闪电哨声波和非闪电哨声波图像特征的区分度，使得分类器识别闪电哨声波的能力增强.

图21 特征分布趋势图Fig.21 Character distribution

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.

4.2.2 特征可分的角度

为了定性地比较不同特征的可分性，首先采用不同的特征提取方式提取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 不同特征训练的分类器的AUCFig.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

5 结论

本文在张衡一号卫星的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接收机观测数据开展闪电哨声波联合观测研究.

致谢本文特别感谢来自应急管理部国家自然灾害防治研究院的张衡一号卫星团队的所有成员为本研究提供的技术服务支持.