2013, Vol. 56
2. 中国极地研究中心, 上海 200136
2. Polar Research Institute of China, Shanghai 200136, China
极光亚暴是太阳风-磁层-电离层耦合系统中最频繁的扰动形式，是地球空间最重要的能量输入、耦合和耗散过程.亚暴的发生可能造成高纬度地区无线电通讯中断，同时导致地球同步轨道附近相对论电子通量的增长[1-2]，从而造成卫星的充电等效应.因此，对亚暴相关问题的研究一直都是日地空间物理学最受重视的核心前沿问题之一.
准确界定亚暴起始时刻是理解亚暴激发机制的一个关键问题，目前已有大量相关研究，比如导致亚暴开始的物理条件/机制是什么以及控制亚暴发生位置及亚暴大小的因素是什么，这些研究已经取得了很大的进展.目前，判定亚暴起始主要有两种途径：一种是基于各类物理参数，采用模式识别的方法自动识别或人工判定亚暴事件起始时刻. Sutcliffe等人[3]提出采用神经网络的方法从Pi2脉动中检测出亚暴起始时刻；Murphy等人[4]基于小波分析理论分析了五个亚暴事件；Tokunaga等人[5]借助奇异值变换从地面磁力计数据中检测出了全亚暴起始时刻；Kataoka等人[6]通过对地磁脉动指数采用Hilbert-Huang变换来研究亚暴起始时刻.然而，由于地面观测台站分布的局限性，由地面观测数据得到的亚暴起始时刻与真实的亚暴起始时刻之间常常存在一定的时延. Meng和Liou[7]曾指出，这些亚暴特征参数中没有一个与极光亚暴有一一对应关系，尤其是在地磁活跃时期.通过这些特征参数求得的亚暴初始时刻可能比真正的极光亚暴提前或延后，这依赖于信号源外面的信号传播速度.
另一种常见的研究模式是基于人眼观测和一定的语义规则，从大量的极光观测图像中直接判断亚暴事件初始时刻. Akasofu[8]首次提出亚暴概念时利用的是局域的全天空极光图像数据.目前，人们普遍认为全域极光图像是目前唯一能区分亚暴空时效应的一种手段[7]，而且是最可信赖的亚暴研究手段[9].已有基于极光图像上的亚暴分析大多是采用案例分析的模式[8-10]，即通过对几组亚暴事件的详细分析来了解亚暴发生期间的空间地磁环境或解释亚暴发生机制.可是，案例分析中少量的数据使得得出的结论在通用性上大打折扣.因此近年来，基于IMAGE卫星和Polar卫星上的紫外成像仪（Ultraviolet Imager，UVI）采集到的全局极光图像，已经出现了一些统计研究工作[11-14].但是，人工确定极光亚暴起始时刻，尤其是基于一系列极光图片来判断是很繁琐的，而且还容易受实验者主观偏差的影响，从而造成结果的不可靠.因此，亟需开发一种自动的和客观的亚暴事件识别技术.特别是由中国空间科学家提出的夸父计划即将实施，其中夸父B卫星将会对极区的全域极光进行长期的连续观测并带来海量的极光图像，此时自动的亚暴事件识别方法需求就更加迫切.已有的文献中还没有涉及从全局UVI极光图像中自动检测极光亚暴的工作. Kim等人[15]从静态极光图像角度分析极光类型，提出了一种基于分层表述的紫外极光图像检索方法，并明确提出“厚极光”可能是亚暴；他们还指出，通过检测极光图像序列中极光卵的厚度变化可以实现对亚暴的检测，但是并没有给出相关实验结果.
本文旨在从大量的Polar卫星UVI图像中自动地检测出亚暴起始时刻.为了实现这一目标，在图像预处理操作之后，我们首先将原始UVI图像通过网格化转换到磁地方时-地磁纬度（MLT-MLAT）直角坐标系下；然后，用带有空间信息的模糊c均值聚类方法[16]提取每帧UVI图像中的亮斑.最后，找出那些出现局部亮斑突然点亮这一现象的图像帧，并考察每一帧后紧邻的一段时间里亮斑在强度、面积以及极向边界的变化情况，判断突然点亮帧是不是亚暴起始所在.通过与Liou在文献[14]中给出的人工标记对比，实验结果证明我们的方法能有效地从UVI图像中自动检测出亚暴的发生时刻.
2 亚暴介绍及相关研究回顾 2.1 亚暴特征描述极光亚暴概念首次由Akasofu[8]通过分析多个基于地面的极光成像仪的数据提出，他把极光亚暴分为两个阶段：膨胀相和恢复相.随着对极光亚暴的进一步认识，Mcpherron[17]根据前人的研究结果得出，在膨胀相之前应该有一个被称为“增长相”的时期.在增长相期间，磁层从太阳风获取能量，为后面的膨胀相期间的能量释放做准备.亚暴增长相时，夜侧极光弧向赤道向运动，但这些赤道向运动的极光弧很难在全局UVI极光图像上进行识别[13]；亚暴膨胀相起始时，夜侧极光在其赤道向一侧突然点亮，并迅速向极向以及东西向膨胀.本文研究的亚暴起始时刻指的是膨胀期起始时极光突然点亮的时刻.
本文的目的在于设计一种能够自动从UVI图像序列中检测亚暴事件初始时刻的方法，所以，首先我们要明确亚暴过程中极光演变的典型特征有哪些.根据前人的描述[12-14]，亚暴过程中极光在形态特征上主要有如下变化：
（1）亮斑突现：午夜时段附近极光卵局部区域突然点亮.亚暴起始往往发生在66°MLAT和23MLT附近[11, 14].
（2）亮斑增强且迅速膨胀：该亮斑向极向和晨昏两侧膨胀（同时伴随轻微的赤道向膨胀），同时亮斑的强度快速增加，直至达到最大.若极光活动较弱，如伪暴，亮斑区域可能不能到达极光卵的极向边界.
（3）膨胀消退：亮斑膨胀到极向最大位置后，亮斑强度和膨胀区域立即开始消退或者过一段时间（几十分钟）开始消退.亮斑膨胀达到的极向最大纬度与亚暴的活动强度有关，强度越大，膨胀的范围越大，甚至可达到80°MLAT以上.
（4）亮斑消失：亚暴结束后，之前出现亮斑的区域的极光强度恢复到亚暴之前的强度水平.以上四种形态变化中，亮斑突现与增强膨胀两过程对应亚暴的膨胀相，膨胀消退与亮斑消失两过程对应亚暴的恢复相.
2.2 相关统计研究文献[11]中，Frey等人通过对IMAGE卫星头两年半航行期间（2000年3月-2002年12月）采集到的UVI极光的考察，共标记了2439个亚暴的起始时刻和位置.这一研究后来扩展到标记了整个IMAGE卫星运行期间的4193个亚暴[12].他们采用的准则如下：一是要求在极光卵中必须能观察到明显的局部极光点亮现象；二是局部点亮必须扩展到极光卵的极向边界并且在磁地方时方向有至少20 min的持续时长；三是两个独立的亚暴之间间隔必须超过30 min以上.其中第二个准则用于消除那些没有演变成亚暴的伪暴，第三个准则消除了相邻亚暴间隔特别近的情况以及多亚暴情况.
根据北半球Polar卫星采集到的UVI图像，Kullen和Karlsson[13]也做了一份统计工作：在1998年12月至1999年2月这三个月的时间里，他们共标记了390个伪暴事件和484个亚暴事件.每一个能在Polar卫星UVI图像上观测到的具有亚暴类似行为特征的事件都考虑进来了，包括那些发生在极光卵晨昏两侧的事件.只有发生在10-14 MLT期间的日侧点亮现象不考虑.他们把那些极光只点亮却没有全局膨胀的事件以及那些有至少一部分时长是在亚暴主要阶段之外的事件定义为伪暴.
Liou[14]对Polar卫星整个在航期间的UVI图像数据做了一份更完备的统计工作：从1996-2000年（卫星运行在北半球期间）的UVI图像库中共标记了2003个亚暴事件，从2007年（卫星运行至南半球）的UVI图像库中标记了536个亚暴事件. Liou采用的准则如下：首先对极光点亮强行加入10 min的时间间隔；其次，为了不丢失一些较弱的亚暴事件，并不要求膨胀区域达到极光卵的极向边界，但是极向膨胀太少（1°-2°MLAT）的事件忽略.最低要求是亮斑在极向方向上的膨胀范围要达到几个地磁纬度，并且向两侧膨胀1~2个磁地方时区域.第三，忽略那些常常在亚暴之后几分钟处发生的极向边界点亮现象.但如果点亮位置是在赤道附近的低纬处，考虑当作一个新的亚暴.如果两个亚暴同时在不同的位置处发生，考虑比较强的那个.
由上面的描述我们可以得知，极光作为亚暴活动的一种表现形式，当亚暴发生时，极光到底会有什么样的特征出现，目前还只有定性的描述，还没有一致公认的量化的界定.
3 数据与方法 3.1 数据集本文使用的极光数据是Polar卫星在1996年12月-1997年2月期间采集的UVI极光图像.采用冬季数据是为了避免日晖效应对亮斑判定造成干扰.紫外成像仪工作时有四个滤波通道，可是，不同的滤波器之间的图像不能直接拿来比较.所以，为了图像的一致性，和文献[13]中的处理方式一样，我们也只采用LBHL（Lyman-Birge-Hopfield long）这个滤波通道下的图像，它的波长在160~180nm之间，图像的空间分辨率约为30 km，每幅图像大小为228行200列.事实上，只有Polar卫星在北半球，并且位于远地点附近的轨道高度，才能看到全域的极光图像.在后面的实验过程中，我们要求UVI图像成像正常，实际上也就是主要考虑的是远地点卫星图像.
3.2 图像预处理在检测亚暴事件之前，和文献[18]中类似，我们先对噪声较强的UVI图像做了如下预处理工作：
（1）移除大范围的背景：根据文献[8]对亚暴的描述，我们把地磁50°以下的极光视为背景，因为亚暴主要发生在高纬区域.移除背景区域的极光有利于重点关注有可能发生亚暴的有效区域.
（2）掩膜：同样为了重点突出有效数据，根据UVI图像外边界呈椭圆状的特点[19]，我们用直径为228像素的圆去截取原图，并从中间把宽度限制到200个像素，最后图像大小依旧为228行200列.
（3）负值点清零：清除因仪器噪声等导致的图像像素值为负值的点.
（4）中值滤波：用于平滑图像.
3.3 MLT-MLAT网格图像转换通过前面极光亚暴发生期间极光形态的分析可知，其主要形态变化集中在亚暴期间凸起亮斑区域的亮度变化及面积变化.与全天空极光图像相比，全局UVI图像呈现出的极光形态特征不是很明显.另外，亚暴期间，凸起主要向极向和晨昏两侧膨胀变化，而这些方向的变化在UVI图像中呈现为一种朝极光卵中心和沿着极光卵外围边界的变化.由于不同UVI图像采集时卫星的位置会发生改变，因此不同UVI图像中同一位置区域并不一定对应着相同的磁地方时和地磁纬度，这就意味着我们不能直接用亮斑在UVI图像间的变化情况来反映其在MLT与MLAT方向的改变.
为了便于从图像中直接识别极光沿极向和晨昏两侧的形态变化，我们对不同时刻采集到的UVI图像进行网格化处理，将指定范围内的磁地方时和地磁纬度扇区按照其MLT和MLAT两个方向铺展为一张平面图，该平面图像称为MLT-MLAT网格图像（x轴表MLT，y轴表MLAT）. 图 1给出了一幅UVI图像及其对应MLT-MLAT网格图像的例子.由于极光亚暴发生在夜侧极光区域，因此，在后续处理中我们采用的MLT-MLAT网格图像均是夜侧时段的图像，如图 1d所示. UVI图像转换为MLT-MLAT网格图像步骤为：
|
图 1 UVI图像与MLT-MLAT网格图像转换示例（1996年12月18日05:36:08） （a）为原始UVI图像；（b）为UVI图像对应的MLT-MLAT网格图像；（c）为（b）平滑处理之后的结果；（d）为（c）的夜侧图像. Fig. 1 An example of transforming UVI images into MLT-MLAT grid images. The image occurred at 05:36:08 of December 18, 1996 (a) The original UVI image; (b) The MLT-MLAT grid image with the x axis being mlt and y axis being mlat; (c) The smooth processing result of (b); and (d) only gives the nightside image. |
（1）设定磁地方时和地磁纬度范围.由于亚暴是一种夜侧极光现象，我们只考虑18-24MLT及00-06MLT范围内的数据；极光主要发生在极光卵的高纬区域（地磁纬度66°附近），因此地磁纬度被限制在50°-90°MLAT；
（2）设定网格大小，Δmlt与Δmlat.本文采用的参数值分别为0.1和0.2.转换后的网格图大小为200×240（（90°MLAT-50°MLAT）/0.2=200行，24小时/0.1=240列），如图 1b；若只考虑夜侧时段则图像大小为200×120（图 1d）.
（3）对网格图里的每个坐标点（mlt，mlat），找到其对应在UVI图像中的具体位置，并计算UVI图像中落在每个网格（mlt：mlt+Δmlt，mlat：mlat+Δmlat）里的像素的均值，该均值即代表了网格图中该网格点的强度.由于Polar卫星UVI成像仪的空间分辨率和成像区域的限制，可能会出现某些网格点在原始UVI图像中对应区域里并没有正常像素值，此时，往四周各方向均扩大一个网格单位区域的查找范围（区域大小由1×1变3×3），再求均值.
（4）平滑滤波.为了便于后续处理，我们对步骤（3）得到的图像进行均值滤波平滑图像.
由MLT-MLAT网格图像的定义知，极光亚暴过程中凸起亮斑区域的形态变化在MLT-MLAT网格图像中主要表现为水平方向和竖直方向的形态变化. UVI图像中的亮斑凸起区域在网格图像中也呈现为一个凸起.各网格图像相同位置对应相同的MLT、MLAT值，因此前后帧间的形态变化可以通过该图像直接表述.值得注意的是，虽然由原始UVI图转成网格图后，极光卵凸起亮斑区域的形状会发生改变；但亚暴过程中我们主要关心亮斑的磁地方时、地磁纬度随时间的变化情况，而把UVI图像转换成网格图像并不影响这些信息的提取.
3.4 基于SFCM算法的亮斑提取聚类分析是数理统计中的一种多元分析方法，它定量地确定样本之间的相似程度/亲疏关系，从而客观地划分样本类型（对样本分类）.考虑到UVI图像中亮斑区域的边界模糊性，本文采用模糊聚类方法来提取极光亮斑.在模糊聚类方法中，模糊c均值（Fuzzy C-means，FCM）[20-21]应用广泛，它通过最小化目标函数（式（1））得到每个样本点对所有类中心的隶属度，从而决定样本点的类属以达到自动对数据样本进行分类的目的：
|
(1) |
式中，n为样本点数，c为聚类数，uij为第j个样本xj相对于第i个聚类中心vi的隶属度，m是控制模糊度的模糊权重，‖*‖是欧式距离表相似度，为了方便下面我们将‖xj-vi‖2简记为di j.对所有输入参量求导，使式（1）达到最小的必要条件为：
|
(2) |
和
|
(3) |
由上述两个必要条件，FCM算法求解是一个简单的迭代过程.当maxij‖uij-ûij‖＜ε时迭代终止，其中常数因子ε∈[0 1].
图像数据的一个很重要的特征是相邻像素之间相关性.换句话说，相邻像素的各类特征值更接近，隶属同一类别的可能性也更大.空间关系在聚类的过程中非常重要，但常规的FCM算法并没有利用空间信息.鉴于此，文献[16]提出了一种包含空间信息的FCM算法（SFCM），其中，空间函数定义为：
|
(4) |
其中NB（xj）指以xj为中心的邻域，本文中采用3×3窗.跟隶属度函数一样，空间函数hij也表示像素xj属于第i个聚类的概率.如果某一像素邻域里的像素都属于同一个聚类中心，这个像素隶属于这一聚类中心的空间函数值就很大.空间函数嵌入到隶属度函数的方式如下：
|
(5) |
其中p和q是控制隶属度函数和空间函数的相对重要性的参数.
聚类有效性（确定聚类数目）是聚类问题中的一个经典问题，在过去几年里被广泛研究[22-23].直观来看，对我们研究的亚暴检测问题，分成三类似乎是最好的，即将UVI图像划分为亮斑区、极光卵区和背景区.但是，通过实验发现，模糊聚类成三类的结果与人眼视觉的判断并不相符.如图 2所示，分成三类得到的亮斑区域（白色区域）并不是人眼判断的亮斑区.在对大量的极光图像进行分析后，我们选择聚类数为6.首先用SFCM算法把每一帧UVI图聚成6类（按强度从强至弱标为类1、类2、…、类6），然后进行区域合并判断出亮斑区域.区域合并的准则如下：
|
图 2 SFCM聚类确定亮斑区域第一列给出的是UVI图像.第二列给出的是SFCM聚成三类的结果，白色、灰色、黑色区域分别表示亮斑区、极光卵区和背景区.第三列给出的是SFCM聚成6类的结果，强度由强至弱分别对应类1（白色区域）到类6（黑色区域）.绿色曲线包围的区域就是本文算法确定的亮斑区域. （a）1997-01-28-12:03:13；（b）1996-12-04-00:02:49；（c）1997-01-04-05:29:58；（d）1997-01-26-22:39:39. Fig. 2 Determining bright bulge areas with SFCM clustering algorithm The left column illustrates the UVI images.The middle column is the clustering results with cluster number c=3, where the white areas, gray areas and black areas represent the bright bulge, polar oval and background respectively.The right column gives the clustering results with cluster number c=6, where the clusters from 1 to 6 are mapped to the clusters with the luminosity from strong to weak and the identified bright bulge areas are labeled with green line. |
令各个聚类中心的强度分别为v1，v2，…，v6.如果v3＞80，则定义亮斑区域包括类1、类2和类3；如果v2＞80，则亮斑区域包括类1和类2.以上两种情况对应于极光强度非常强的情形.对于一般的极光强度，我们找出相邻聚类中心值的差值中的最大值，如果v1-v2值最大，则定义亮斑区域只包括类1；如果v2-v3值最大，定义亮斑区域包括类1和类2；如果v3-v4值最大，定义亮斑区域包括类1、类2和类3.如果上述方法得到的亮斑区域包含多个独立的块，我们只考虑面积最大/强度最强的那个块（图 2d所示）.从图 2给出的例子看，经过区域合并得到的亮斑区域和人眼的判断基本一致.
3.5 亚暴初始时刻确定已有的统计工作大都是每位作者在自己制定的互不相同的准则下完成的.由于本文考虑的时间段只有Liou[14]的标记里包含到了，为了更方便地评价算法的性能（通过将检测结果与人工标记做对比），本文采用文献[14]中提出的那些准则.具体步骤如下：
（1）检测哪些帧中包含了突然点亮的局部极光（凸起）.首先要求图像成像必须正常，忽略那些在亚暴发生有效区域里是全黑图像或者成像很不全的情况（因为主要考虑亮斑在极向和晨昏两侧的膨胀，因此“成像不全”指的是亮斑整个外围或亮斑在上、左、右三个方向的任意一个方向上有超过三分之一的地方是紧挨着全黑区域）.既然是亮斑，我们要求极光强度不能太低（大于25）.关于“突然点亮”，我们指的是亮斑区域极光强度不仅比周边区域强，而且也比前一帧同一位置处的极光强度强（强度值绝对变化大于15或相对变化大于1.4倍）.
（2）在每一幅包含局部极光突然点亮的UVI图像之后，我们继续考察其后20 min左右的极光图像序列，综合判断当前事件是不是一个亚暴.在考察的这段图像序列中，如果亮斑消失或者亮斑位置在相邻帧间发生了大的改变（前者意味着当前亚暴膨胀期已结束，后者意味着当前帧已经属于下一个亚暴事件了），考察序列时长缩短至出现这些现象的前一帧图像结束.相反，若第20 min处的图像还处在亮斑膨胀期的关键时刻，则适当延长考察序列时长（为了避免进入下一个亚暴时间而造成判断失误，最长时长限制在30 min以内）.
（3）对上一步骤挑选出来的极光突然点亮后的待考察极光图像序列段（~10-30 min），我们主要通过考虑亮斑在强度、面积及极向边界的变化情况来判断当前事件是不是一个亚暴事件.由于亚暴膨胀期亮斑的改变是波动式而非单调连续的[24-25]，实验时我们记录考察序列段中亮斑所达到的最大强度、最大面积、极向边界最高纬度值，并与点亮帧中初始亮斑的强度、面积、极向边界比较.如果增幅较大，我们判断是亚暴事件，点亮帧所在的时刻就是亚暴初始时刻；反之如果增幅很小，则忽略不计.具体实验过程中，我们要求极向边界往极向方向移动超过1.2°MLAT，膨胀面积超过20个坐标像素，亮斑强度增幅超过8.由于有些亚暴过程中的关键极值时刻可能没有被LBHL滤波器采样到，所以，在本文数据集上来进行亚暴识别存在一定的局限性.从我们所用的图像数据库中看，很多Liou标记的亚暴事件不能同时满足上述亚暴条件.经过分析，我们针对一些常见的情况放松了条件：当强度增强较大（大于40）且面积增大较多（大于35）时，极向移动可以较少（大于0.6°MLAT）；或者如果亮斑是由很暗突变到很亮（点亮帧亮斑强度超过前一帧的3倍或强度值增幅超过30），则允许亮斑强度在随后图像序列中不再增加.
（4）在出现突然点亮的亮斑之后，可能会多次出现亮斑强度变弱后又再次增强或亮斑面积缩小后又再次极向膨胀的情况.出现这种情况时，我们只考虑第一个，忽略那些位于膨胀期的极光再点亮和再次膨胀的情况.参考Liou的标记，我们要求两个亚暴初始之间间隔超过12 min.
4 实验结果及分析 4.1 数值评估实验共检测出489个亚暴事件，其中有242个与Liou的标记一致（在这期间，Liou共标记了262个亚暴事件）.当然，在这些一致的结果里，可能会有几分钟的出入.这是因为，一方面，由于UVI图像时间分辨率不高，相同波段UVI图像的最小采样间隔为半分钟，一般是~3 min，很难精确界定亚暴起始时刻[13, 26]；另一方面，Liou的亚暴初始点判断是采用了多个滤波通道数据得到的.实验结果详见表 1.
|
表 1 亚暴事件检测结果 Table 1 Experimental results of the substorm event |
从表 1可以看出，这三个月中亚暴查全率均高于90%，也就是说Liou标记的亚暴事件中绝大多数都被检测出来了；可与查全率相对应的查准率只有50%左右，也就是说，我们检测的事件中有差不多一半Liou没有标记.深入分析Liou的标记和本文实验结果，我们发现Liou的标记很不全，而且还受人为偏差的影响.更详细的结果分析见4.2节.
4.2 结果分析本文方法检测到的亚暴结果和Liou的标记之间存在一些差异.不仅一些检测到的亚暴起始点与标记的相比有几分钟的出入，还存在少量的标记事件被漏检和大量的多检事件（即Liou未标记但本文方法检测到了的事件）.
4.2.1 漏检亚暴分析图 3列举了一些Liou标记为亚暴，但本文方法漏检了的亚暴事件例子.其中，第一个例子（a）漏检的原因是强度太弱，视觉上已经看不出有亮斑的存在；第二个例子（b）漏检的原因是突亮强度弱，而且极向膨胀太少；第三个例子（c）漏检的原因是相邻图像之间间隔太长造成快速演变的亚暴事件关键细节丢失，亚暴膨胀期只采集到了强度最强的那一帧.
|
图 3 漏检亚暴事件举例 （a）1996-12-21-02:32:07；（b）1997-01-02-18:42:38；（c）1997-02-14-01:04:17. Fig. 3 Some examples of missed substorm events |
当确定亚暴初始点时，Liou参考了LBHL和LBHS两个波段的紫外全天空图像并分别得到一个初始时刻点，然后综合考虑确定亚暴初始为那两个时刻中较早的那个或是二者之间的某一时刻.但是本文实验中只使用了LBHL这一个波段的数据.这一方面导致我们的实验结果与标记的时刻有些小差异（常常是略微延迟），因为亮斑初始点亮的瞬间可能没有被LBHL滤波器采样到；另外更严重地是，本文所用数据集中相邻帧之间的时间间隔基本都在~3 min，这么长的间隔会漏掉很多有用信息，特别对那些极光演变特别迅速的亚暴事件来说，常常出现膨胀期从开始至结束这一过程只捕捉到了一帧图像或者没有捕捉到膨胀期的最强亮斑/最大面积/极向最高边界所在的时刻（如图 3c所示），此时，就不可避免地造成了这些亚暴事件的漏检.分析所有漏检的21个亚暴事件，我们发现原因都基本类似：即从LBHL波段的图像序列中我们看不到突然点亮的亮斑强度增强和极向膨胀这些亚暴发生时必然会出现的最基本特征，也就是说从图像序列的角度我们不能判断当前事件是一个亚暴.
4.2.2 多检事件分析本文实验检测到结果是那些UVI图像序列中亮斑变化具有亚暴特征的极光事件，仔细分析这些Liou没有标记的多检事件，我们可以将其大致分为如下四种情况.
（1）是亚暴事件而Liou漏标了. 图 4给出了发生在1997年1月的四个例子.从图 4给出的UVI图像序列，我们很容易看出这些事件均满足亚暴定义中对亮斑变化的那些要求；此外，我们还发现这四个事件也被包含在从SuperMAG求得SME指数确定的亚暴事件列表中[27-28].
|
图 4 多检事件举例之是亚暴事件 （a）1997-01-03-21:49:58；（b）1997-01-11-11:40:33；（c）1997-01-28-20:33:31；（d）1997-01-29-05:25:16. Fig. 4 Examples of some extra identified events which are substorm events |
（2）属于膨胀期再点亮过程.已有文献[24-25]指出，亚暴膨胀期的变化通常是步进式的而非单调连续的，文献[25]还通过一个发生在1997年5月16日的亚暴事件中亮斑中心点的位置变化来说明这一点.这导致人们常常很难，尤其是算法很难判断极光重新点亮的过程是属于前一亚暴的膨胀期还是新的亚暴事件.本文多检结果中包含了少数需要被剔除的属于膨胀期的事件.事实上，在Liou的标记中，我们也发现有些事件可能是前一亚暴的膨胀期极光再点亮.比如在1996年12月11日这一天Liou标记的发生在09:54:36和10:07:22两事件，从UVI图像序列中看二者更可能就是同一个亚暴事件.这段时间里由SME指数得到的亚暴列表[27-28]也只包含了一个亚暴事件.类似地，1997年2月20日这一天Liou标记的发生在20:10:39和20:26:26的两事件，从UVI图像序列里看也更像是同属一个亚暴事件.Frey[11-12]在标记亚暴的时候，是通过要求两相邻亚暴时间间隔必须在30 min以上来减少对这种情况的判断.
（3）伪暴情况.目前为止，还没有定量的标准来区分亚暴、极向边界增强（PBIs）和伪暴这三种极光现象.文献[29-30]已经分别指出亚暴不能和伪暴、PBIs区分开.一部分研究者[13, 31-32]认为伪暴就是弱的亚暴事件.文献[33]也指出一般情况下很难将伪暴（即使是孤立的或增长期的伪暴）、亚暴以及PBIs区分开，因为这三类极光增强事件呈现出相似的极光特征和地磁特征.在Peter D. Boakes的博士论文[34]中，他指出Frey的亚暴标记[11-12]中可能包含了一些并非真正的亚暴事件，而且还可能包含了大量更弱的事件，这些事件可能会带来和亚暴相似的极光点亮，如伪暴、PBIs等.Liou的亚暴定义标准要比Frey要求的条件宽松，因此Liou的标记中也肯定包含了一些伪暴事件.本文实验中的参数是根据Liou标记的亚暴事件得出的，所以实验结果中不可避免地会出现一些伪暴事件. 图 5列举了实验结果中三个可能是伪亚暴的事件，这三个事件Liou的标记中没有，可是有一些与之非常类似的情况他又标记了.
|
图 5 多检案例分析之伪亚暴事件 （a）1997-01-01-22:00:11；（b）1997-01-11-13:30:57；（c）1997-01-28-11:17:13. Fig. 5 Examples of some extra identified events which are pseudobreakups |
（4）模棱两可的情况.实验结果中还有个别事件，从图像上很难界定它们是不是亚暴事件，我们称这些为模棱两可的事件，如图 6所示.出现这种情况一方面是因为实验用到的LBHL这个波段数据的采样间隔太长，不能全面准确地反映极光亮斑的变化情况；另外，也是因为亚暴现象在UV I图像上的界定比较模糊所致.
|
图 6 多检案例分析之模棱两可的事件 （a）1997-01-01-01:10:12；（b）1997-01-29-16:25:50. Fig. 6 Examples of some extra identified events which are ambiguous cases |
仔细查看Liou的标记，我们发现有一些非常类似的事件有时被标记成了亚暴事件，有时却又没有被标记.这其实就是当人们通过眼睛来检测极光亚暴时，不同人持有的标准会有所差异，而且由于无法定量化所以即使是同一个人的判断标准也难以始终保持一致，这些情况一方面会造成对那些本身就模棱两可的事件标记的不一致；另一方面，本文方法的参数是通过参考Liou的标记进行调整确定的，这些不一致的标记会干扰参数的设定进而影响实验结果.
5 结论及讨论针对UVI图像中亚暴事件膨胀期初始时刻的判定问题，本文提出了一种基于模式识别的自动检测方法，并在1996年12月-1997年2月这三个月的Polar卫星紫外图像数据上进行了实验验证.实验结果显示，本文的方法能有效地找出具有亚暴特征的极光事件.虽然实验结果与Liou的标记相比查准率不高，但这很大一部分原因是目前极光亚暴特征尚未有明确的界定，同时Liou的标记也存在一些问题等造成的.所以，不能单用查全率查准率来评价本文方法.但不管怎样，高达92%的查全率（8%的漏检率）证明本文的方法还是很有实用价值的：可以在海量的紫外极光图像数据中，有效地完成亚暴事件的初步筛选：这就能大大缩小原始数据集，从而有利于研究人员对极光亚暴现象进行深入的分析与探讨.与人工标记方法相比，本文的方法具有标准统一、省时省力等优点.
可是，由于亚暴过程非常复杂，单纯从极光图像来完成自动识别存在一定难度.对本文方法初步筛选后得到的候选亚暴事件，我们还需人工去一一识别到底是不是亚暴.特别是对于亚暴、极向边界增强（PBIs）和伪暴这三种比较相似的极光现象更需要人工进行辨识.在分析亚暴发生发展的细节特征时，我们除了监测极光的活动特征外，还需要分析AE指数的变化、Pi2的活动情况、以及磁尾磁场的偶极化和相应的亚暴能量粒子注入情况等.
极光亚暴初始时，UVI图像上表现出来的特征目前还只有定性的描述，尚未有公认的定量的界定.本文方法在做量化处理时涉及很多参数值，如什么情况算亮斑出现，亚暴现象要求亮斑至少需要极向膨胀多大范围及强度增加多少，本文中这些参数是基于1996年12月-1997年2月这三个月里Liou标记的亚暴事件，权衡考虑查全率和查准率后得到的.当数据集/人工标记发生了改变或者若实际应用中对漏检/多检有明确的偏好要求时，这些参数值需要重新调整.另外，由于本文所用数据的时间分辨率低，很多有用信息被丢失，部分亚暴事件在我们用到的数据集上不能同时满足亚暴定义中的亮斑强度增加和极向膨胀这些条件.因此，在分析了大量亚暴事件后，实验中我们放松了亚暴条件，比如，当突然出现极光亮斑强度特别强时（说明亮斑点亮的真实时间介于当前帧和前一帧两个时刻之间），我们并不要求这个亮斑强度再增强（因为当前帧可能已经是最强亮斑所在了）.如果上述这两个问题解决了，即当数据集时间分辨率足够高而且极光亚暴初始特征有明确的界定时，我们的方法就能更为准确地找出亚暴膨胀期的初始时刻.
后续工作中，针对Liou标记存在一些误差的问题，我们将尝试换其它如Frey或Kullen等人的标记，进一步验证本文方法的性能. UVI成像仪在1996年冬季是多个滤波器同时工作的，所以当我们只用LBHL这一个波段数据时就会出现因分辨率低而造成决策困难的情况.下一步我们将尝试分辨率高一些的其它时段的极光数据，以获取更准确的结果.
致谢非常感谢美国国家空间科学和技术中心（The National Space Science and Technology Center，NSSTC，U.S.）提供的UVI图像数据以及Liou提供亚暴事件的人工标记.
| [1] | 李柳元, 曹晋滨, 周国成. 磁层相对论电子通量变化与磁暴/亚暴的关系. 地球物理学报 , 2006, 49(1): 9–15. Li L Y, Cao J B, Zhou G C. Relation between the variation of geomagnetospheric relativistic electron flux and storm/substorm. Chinese J. Geophys. (in Chinese) , 2006, 49(1): 9-15. |
| [2] | Li L Y, Cao J B, Zhou G C, et al. Statistical roles of storms and substorms in changing the entire outer zone relativistic electron population. J. Geophys. Res. , 2009, 114(A12): A12214. DOI:10.1029/2009JA014333 |
| [3] | Sutcliffe P R. Substorm onset identification using neural networks and Pi2 pulsations. Ann. Geophys. , 1997, 15(10): 1257-1264. DOI:10.1007/s00585-997-1257-x |
| [4] | Murphy K R, Rae I J, Mann I R, et al. Wavelet-based ULF wave diagnosis of substorm expansion phase onset. J. Geophys. Res. , 2009, 114: A00C16. DOI:10.1029/2008JA013548 |
| [5] | Tokunaga T, Yumoto K, Uozumi T, et al. Identification of full-substorm onset from ground-magnetometer data by singular value transformation. Mem. Fac. Sci., Kyushu Univ., Ser. D, Earth & Planet. Sci. , 2011, 32(3): 63-73. |
| [6] | Kataoka R, Miyoshi Y, Morioka A. Hilbert-Huang Transform of geomagnetic pulsations at auroral expansion onset. J. Geophys. Res. , 2009, 114(A9): A09202. DOI:10.1029/2009JA014214 |
| [7] | Meng C I, Liou K. Substorm timings and timescales: A new aspect. Space Sci. Rev. , 2004, 113(1-2): 41-75. DOI:10.1023/B:SPAC.0000042939.88548.68 |
| [8] | Akasofu S I. The development of the auroral substorm. Planet. Space Sci. , 1964, 12(4): 273-282. DOI:10.1016/0032-0633(64)90151-5 |
| [9] | Liou K, Meng C I, Lui A T Y, et al. On relative timing in substorm onset signatures. J. Geophys. Res. , 1999, 104(A10): 22807-22817. DOI:10.1029/1999JA900206 |
| [10] | Liu J M, Zhang B C, Kamide Y, et al. Spontaneous and trigger-associated substorms compared: Electrodynamic parameters in the polar ionosphere. J. Geophys. Res. , 2011, 116(A1): A01207. DOI:10.1029/2010JA015773 |
| [11] | Frey H U, Mende S B, Angelopoulos V, et al. Substorm onset observations by IMAGE-FUV. J. Geophys. Res. , 2004, 109. DOI:10.1029/2004JA010607 |
| [12] | Frey H U, Mende S B. Substorm onsets as observed by IMAGE-FUV.//Proceedings of Eighth International Substorm Conference, Univ. of Calgary, Banff Centre, 2006: 71-76. |
| [13] | Kullen A, Karlsson T. On the relation between solar wind, pseudobreakups, and substorms. J. Geophys. Res. , 2004, 109: A12218. DOI:10.1029/2004JA010488 |
| [14] | Liou K. Polar ultraviolet imager observation of auroral breakup. J. Geophys. Res. , 2010, 115(A12): A12219. DOI:10.1029/2010JA015578 |
| [15] | Kim S K, Ranganath H S. Content-based retrieval of aurora images based on the Hierarchical Representation. Advanced Concepts for Intelligent Vision Systems, Part II, LNCS , 2010, 6475: 249-260. DOI:10.1007/978-3-642-17691-3 |
| [16] | Chuang K S, Tzeng H L, Chen S, et al. Fuzzy c-means clustering with spatial information for image segmentation. Computerized Medical Imaging and Graphics , 2006, 30(1): 9-15. DOI:10.1016/j.compmedimag.2005.10.001 |
| [17] | Mcpherron R L. Substorm related changes in the geomagnetic tail: The growth phase. Planet. Space Sci. , 1972, 20(9): 1521-1539. DOI:10.1016/0032-0633(72)90054-2 |
| [18] | Wang Q, Meng Q H, Hu Z J, et al. Extraction of auroral oval boundaries from UVI images: A new FLICM clustering-based method and its evaluation. Advances in Polar Science , 2011, 22(3): 184-191. |
| [19] | Cao C G, Newman T S, Germany G A. New shape-based auroral oval segmentation driven by LLS-RHT. Pattern Recognition , 2009, 42(5): 607-618. DOI:10.1016/j.patcog.2008.08.018 |
| [20] | Bezdek J C. Pattern Recognition with Fuzzy Objective Function Algorithms. New York: Plenum Press, 1981 . |
| [21] | Bezdek J C, Ehrlich R, Full W. FCM: the fuzzy c-means clustering algorithm. Computers & Geosciences , 1984, 10(2-3): 191-203. |
| [22] | Sun H J, Wang S G, Jiang Q S. FCM-based model selection algorithms for determining the number of clusters. Pattern Recognition , 2004, 37(10): 2027-2037. DOI:10.1016/j.patcog.2004.03.012 |
| [23] | Wang Y N, Li C S, Zuo Y. A selection model for optimal fuzzy clustering algorithm and number of clusters based on competitive comprehensive fuzzy evaluation. IEEE Transactions on Fuzzy Systems , 2009, 17(3): 568-577. DOI:10.1109/TFUZZ.2008.928601 |
| [24] | Kisabeth J L, Rostoker G. The expansive phase of magnetospheric substorms: 1. Development of the auroral electrojets and auroral arc configuration during a substorm. J. Geophys. Res. , 1974, 79(7): 972-984. DOI:10.1029/JA079i007p00972 |
| [25] | Liou K, Meng C I, Wu C C. On the interplanetary magnetic field By control of substorm bulge expansion. J. Geophys. Res. , 2006, 111(A9): A09312. DOI:10.1029/2005JA011556 |
| [26] | Voronkov I O, Donovan E F, Samson J C. Observations of the phases of the substorm. J. Geophys. Res. , 2003, 108(A2): 1073. DOI:10.1029/2002JA009314 |
| [27] | Newell P T, Gjerloev J W. Evaluation of SuperMAG auroral electrojet indices as indicators of substorms and auroral power. J. Geophys. Res. , 2011, 116(A12): A12211. DOI:10.1029/2011JA016779 |
| [28] | Newell P T, Gjerloev J W. Substorm and magnetosphere characteristic scales inferred from the SuperMAG auroral electrojet indices. J. Geophys. Res. , 2011, 116(A12): A12232. DOI:10.1029/2011JA016936 |
| [29] | Ohtani S, Anderson B J, Sibeck D G, et al. A multisatellite study of a pseudo-substorm onset in the near-earth magnetotail. J. Geophys. Res. , 1993, 98(A11): 19355-19367. DOI:10.1029/93JA01421 |
| [30] | Rostoker G. On the place of the pseudo-breakup in a magnetospheric substorm. Geophys. Res. Lett. , 1998, 25(2): 217-220. DOI:10.1029/97GL03583 |
| [31] | Nakamura M, Baker D N, Yamamoto T, et al. Particle and field signatures during pseudobreakup and major expansion onset. J. Geophys. Res. , 1994, 99(A1): 207. DOI:10.1029/93JA02207 |
| [32] | Aikio A T, Sergeev V A, Shukhtina M A, et al. Characteristics of pseudobreakups and substorms observed in the ionosphere, at the geosynchronous orbit, and in the midtail. J. Geophys. Res. , 1999, 104(A6): 12263. DOI:10.1029/1999JA900118 |
| [33] | Kullen A, Karlsson T, Cumnock J A, et al. Occurrence and properties of substorms associated with pseudobreakups. J. Geophys. Res. , 2010, 115(A12): A12310. DOI:10.1029/2010JA015866 |
| [34] | Peter D B. Investigating the relationship between open magnetic flux and the substorm cycle[Ph.D.thesis]. Radio and Space Plasma Physics Group, University of Leicester, 2009. |