极光是人们能够用肉眼直接观测到的具有极区特征的高空大气物理现象，是极区日地物理过程(特别是磁层-电离层相互作用)最集中的表现形式.目前，极光数据的获取有多种方式，包括极光雷达观测、全天空成像、卫星成像、地磁观测、二维宇宙噪声成像观测以及子午面扫描光度计观测等.本文研究利用全天空图像序列，实现对极光活动的连续观测和变化的监视.

在海量数据中准确、快速、自动地检测极光事件是进一步研究极光物理机制的基础，是有效利用日积月累产生的海量极光数据面临的首要问题.为了检测极光活动的变化，快速浏览极光活动的变化情况，空间物理研究人员通常采用的一个方法是制作Keogram图(Yang et al.，2000).它是从一段的全天空极光视频的每一帧极光图像中抽取沿地磁南北方向一条极光亮度信息，以时间为序列排列成极光活动图.显然，Keogram图方法丢失了极光二维的形态和运动的大量信息，不能完整地反映复杂多变的极光活动.通过Keogram图浏览极光活动，会错过许多重要的极光形态和运动信息，例如极光弧上小尺度涡旋结构，以及发生在全天空图像天顶南北磁子午线以外的极光事件等等.

现阶段，人们通过人工手动和计算机自动分析两种方式展开对极光形态和动态特征的研究.在人工分析方面，研究人员通过手工标记和摘选，对极光现象进行逐一研究和分析(Haerendel et al.，1993；Kimball et al.，1998； Hu et al.，1999；Blixt et al.，2005；Qiu et al.，2013).例如，前人往往是通过人眼观察跟踪极光的某些特征，进而估计极光的运动速度和方向.传统的人工研究方法存在很多不足之处.首先，这种方法加入了过多的人为主观认识.例如极光“应该”提取的特征点是由研究人员主观选择的、算法的评估也往往是由人主观决定的；其次，手工标记的特征是通过人工主观筛选得到的元素级别的特征，这些特征缺乏对宏观特征的描述，而且对成像设备的依赖程度很高；最后，由于人工分析的效率不高，得到研究的极光事件的数量大大受到了限制，大多数对极光的研究分析被局限于案例分析(Case study)，在案例分析中只有非常少数的几个短暂的极光事件能够得到研究(Marklund，1984，Kauristie et al.，2001).综上所述，传统的方法一方面对于浩瀚的极光数据来说只能是杯水车薪，另一方面使得到的结果难以得到推广.

在借助计算机自动分析方面，现阶段无论是基于单帧图像的静态分析或者是基于序列的动态分析，所有的极光图像或者视频的表征方法都是基于整幅极光图像的(Wang et al.，2010；Yang et al.，2011；Yang et al.，2012a； Yang et al.，2012b； Yang et al.，2013).这就会带来两个问题.首先，由于同一时刻天空中可能会发生多个极光事件，那么出现在同一幅极光图像中的极光可能属于不同的类型，由不同的磁层过程引起，那么区分图像中的对象和极光图像是有必要的.其次，人们越来越关注小尺度的极光结构，例如褶皱和涡旋，它们具有非常重要的研究价值.然而，现有的基于宏观特征的自动分析方法根本无法处理这些小尺度极光结构.

本文的目标在于量化极光的形态和运动的变化情况，开发一种检测极光活动变化的方法，用于快速监视极光活动.所以，我们需要一种方法能够全面的提取极光的形态和运动信息，包括极光形态分布、运动方向和速度等.这些二维信息的变化都是极光活动的重要特征，能够反映以上信息的有效方法是提取视野内的运动向量场.在全天空视野中，极光的全局运动可以用运动向量场描述，它反映了图像上每一像素点的变化趋势，可看成图像平面上运动而产生的瞬时速度场.通过运动向量场可以了解极光活动许多重要的特征，如极光的形状、运动轨迹和速度等.

如上所述，为了能够有效地表征极光的形态和运动，首先要解决的问题是准确提取极光的运动场.Horn等(1981)提出了光流场方法，它是运动场的二维平面模拟.Horn等将光流场求解问题分为两个部分：数据约束方程和正则化方法，这一研究思路奠定了光流算法发展的基石.针对流体的特性，Corpetti等(2002)做出了以下假设：图像的亮度与流体的密度成正比例关系，由此引入了流体动力学中的连续方程作为构成数据约束的基础，结合二阶散度-旋度正则化方法提取出了大气图像的运动场.Blixt等(2006)首次将光流场的方法用于估计极光的运动场，虽然实验结果体现出了光流场方法应用于极光数据的有效性，但是光流场的基本假设不适用于极光数据，这一点在他们的研究中有所提到并且给出了图像中不符合假设的区域.光流场方法的数据约束方程是基于亮度不变假设的，即观测对象的亮度在运动过程当中保持不变.显然，极光在产生和消亡的过程中，不能满足亮度不变假设.根本原因是在经典的计算机视觉研究背景之下，运动物体往往被假设是刚性的，存在稳定突出的轮廓特征，而在大多数极光图像中并不存在这样明显的观测“目标”，而且极光的形状、亮度、体积等在演变过程中会发生变化.在过去的三十多年的发展中，人们提出了很多种运动场的提取方法.然而，已有的运动场估计方法用于极光视频图像主要存在以下两个问题：首先，在极光运动的过程中，引起发光的高能微粒流的体积、密度以及形状都会发生变化.这些粒子流是可压缩的，透明的，或是可能还会消失.所以，传统的基于亮度相关性的运动场估计方法不适用于极光数据；其次，采用一阶正则化方法的运动估计算法可以得到较为平滑的运动场，而采用二阶正则化方法的运动场估计算法更适合旋转运动.现阶段，几乎所有的运动场估计方法都采用多尺度方法解决运动场的非线性估计问题.但是这些研究针对不同尺度的运动采用了完全一样的正则化方法，这就默认极光在不同尺度下具有相同的运动特征.但是极光的运动极为复杂，不同尺度的运动形式具有不同的特点.所以，估计极光运动采用统一的正则化方法是不合适的.需要指出的是，Verri和Poggio(1987)指出，对于强纹理图像，运动场和光流场是近似相等的.在本文中，我们称估计得到的二维向量场为运动场.

本文旨在基于全天空极光图像的形态和运动情况检测极光活动的变化，为了实现这一目标，我们将构造一条能够检测极光活动变化的曲线，这一曲线的变化基于复杂的极光形态和运动特征，曲线的局部峰值的出现指示极光活动发生变化的时间，从而实现快速浏览极光活动变化的目的.为了反映极光的形态和运动特征，我们首先提取全天空极光视频序列的运动场，然后基于运动场对极光序列进行表征，最后计算滑窗内连续的前后两段视频序列的差异，从而实现量化极光的形态和运动的变化情况，滑动滑窗得到检测极光活动变化的曲线.

本文的观测数据来自我国北极黄河站(地理坐标为(78.92°N，11.93°E)，修正磁纬为76.24°N，磁地方时(MLT)=世界时(UT)+3.1 h)的全天空极光图像序列.为了方便后续操作，我们对原始的采集数据进行了以下预处理操作：减去暗电流、图像灰度拉伸、旋转已经剪裁，详见前期研究工作(Wang et al.，2010).

极光图像中的亮度不同于传统的模式识别研究中使用的图像数据，普通图像中的亮度来自反射光，研究中往往假设目标的亮度不变，极光是高能带电粒子沿高纬磁力线运动，与地球大气层碰撞后产生的地球物理现象，极光的亮度在运动过程中会发生形变、增长以及消亡等.由于极光自身的这些特性，传统的基于亮度不变假设的光流场方法不能应用于极光运动场的估计.

事实上，极光是由太阳风和地球磁场相互作用产生的，期间产生并释放出高达10⁶ MW的能量，这些能量通过电流沿着磁力线流进地球高层大气.太阳风是等离子体，具有磁流体的特性，这些高能带电粒子与地球大气相互作用释放的光子数通量越高，图像的亮度越高.换句话说，单位像素的亮度与极光沉降粒子的能通量密度呈正比例关系.因此，估计极光运动时，本文采用Corpetti等(2002)提出的将流体力学中的连续方程用于运动建模的方法.连续方程的理论基础是质量和动量守恒定律，并且广泛地应用于水波、大气、烟等运动图像的分析.全天空图像理论上提供了地平线以上天空中的所有信息，假设在天顶角80°以内区域的极光亮度与全天空图像的亮度成正比例关系.那么，在基于连续方程构造亮度约束时，即可用极光的亮度代替连续方程中的密度ρ.

其次，考虑到极光运动复杂多变的特点，本文根据北极黄河站观测数据的时间和空间分辨率，将全天空极光图像进行4层塔形分解，记为0~3层(0层为原始图像).根据不同尺度的极光运动特征，采用了多尺度分析方法.利用高斯塔形分解对原始图像进行三层下采样.原始图像(记为level 0)的大小为440×440个像素，Level1到Level3的图像大小依次为220×220、110×110、55×55.这种方式被看作是一种塔形时空运动检测器，其中在粗尺度上检测较大尺度的极光运动，细尺度上检测较小尺度的极光运动.

然后根据观测对象的运动特性决定每层采用的正则化因子，从而使不同尺度下不同的极光运动都得到准确的恢复.对第0、2和3塔形分解层，选择一阶平滑正则化因子；对第1塔形分解层，选择二阶散度-旋度正则化因子.在每一个塔形分解层上，执行基于连续方程的数据约束结合非一致正则化因子的运动场估计方法，全局目标方程为：

得到待检测全天空极光视频序列对应的运动向量场序列V={v(x，t)，x∈Ω，t=1，2，…，T-1}，其中ξ₁是基于连续方程的数据约束方程，ξ¹2是一阶平滑正则化因子，ξ²₂是二阶散度-旋度正则化因子，E(x，t)表示图像的亮度，x指示空间位置，t表示时间，d为需要估计的运动场，以及f₁和f₂分别为二次的惩罚方程.根据以上方法得到全天空极光序列的瞬时运动场，如图 1所示.图 1中的红色箭头的位置表示了极光出现的位置，箭头的方向表示极光运动的方向，箭头的长度描述极光的运动速度大小.因为显示的全天空极光图像的大小为440×440个像素，无法清晰显示每一个像素的运动矢量.为了清晰显示出密集运动场的箭头大小和方向，图 1a是稀疏运动场，即每一个4×4的邻域内显示一个运动矢量.图 1b是1a图蓝色方框中的局部邻域的密集运动场.

图 1

Fig. 1

图 1 全天空极光序列的全局稀疏运动场(a)和局部密集运动场(b)Fig. 1 The global dense motion field and the local sparse motion field of ASI sequence

量化极光活动的变化的第二个关键问题是对极光序列的表征.由于极光形态和运动模式异常多变，在特征空间的表现是类内方差大.对于静态的纹理图片，灰度共生矩阵的方法得到了很好的效果，它是对邻域内像素的灰度值与中心像素的灰度值的共生情况进行统计.受到了灰度共生矩阵的启发，我们将对邻域内某一向量与中心向量的差进行二维直方图统计，统计范围为全天空镜头视野及视频持续时间的时空域内，由此得到一组直方图向量，将不同方向上得到的直方图向量级联起来构成了描述全天空极光视频序列的特征.

首先计算在一种邻域关系下的局部向量差.这里的邻域关系表示在邻域N_x中，选择和中心像素相互距离为d的另一个像素，他们组成了一种邻域关系，d是一个矢量，在这个邻域内，像素与d一一对应，那么就用d来表示这一邻域关系.如图 2所示，N_x选取为在t帧内以x为中心像素的方形邻域，中心像素以及领域内的与中心像素夹角为0的像素(如图 2标记为0°的像素)为一种邻域关系.计算当前邻域关系下的空间邻域N_x内的局部向量差：

图 2

Fig. 2

图 2 基于运动场对全天空极光视频序列的表征方法Fig. 2 Schematic of building the representation of ASI sequence based on motion field

其中x，x+d∈N_x，t≤T-1.统计时空域{Ω×(1，…，T-1)}上所有局部向量差的长度和方向，即可得到当前邻域关系d下的局部向量差的二维极坐标直方图，其中Ω表示全天空极光图像的圆形掩模区域，T表示极光视频序列的长度，x表示像素的空间位置，t表示时间.将在极坐标平面上统计得到的二维直方图记为h^d.

变换不同邻域关系，按照以上方法将会得到不同邻域关系下的时空统计二维直方图.将得到时空统计二维直方图进行级联[h^d₁(t)，h^d₂(t)，…，h^d_m(t)]，其中m是在邻域N_x中选取的邻域关系的个数，d_m表示第m种邻域关系，最后进行归一化操作，如图 2所示.

将一段待检测的极光视频序列S_a输入以上提出的运动场估计器中，首先对每一帧图像建立0~3层塔形分解层.在各分解层上，利用基于连续方程的数据约束，并选择合适的正则化因子，求得运动向量场序列V={v(x，t)，x∈Ω，t=1，2，…，T-1}.

为了检测极光活动的变化，将时间滑窗从视频起始位置滑向结束位置，每次滑动的步长为一帧.取窗宽为2t_w，将第t(t_w< t≤T-t_w，t_w«T)次滑窗内的极光运动向量场序列分为前t_w帧和后t_w帧子序列，得到第t次滑窗内的前半段运动场向量序列V_pre(t)(即在时空域{Ω×(t-t_w-1，…，t)}的运动向量场)的局部向量差时空统计特征，记为f_pre(t).同理得到后半段运动场向量序列V_post(t)(即在时空域{Ω×(t+1，…，t+t_w)}的运动向量场)的局部向量差时空统计特征，记为f_post(t).度量当前滑窗内的前后极光视频序列的差异，计算f_pre(t)和f_post(t)的-square统计距离，将值赋给变化曲线在t时刻的值:

-square统计距离定义为：

图 3

Fig. 3

图 3 用于检测极光活动的变化曲线生成方法的流程Fig. 3 The flow charts of obtaining the change curve for detecting the auroral activities

利用上述方法对2003—2004年越冬观测的极光活动变化进行检测，在以下的实验中我们将输出的变化曲线与对应的Keogram图进行对比.Keogram图是对每一段的全天空极光视频，抽取每一帧极光图像中沿地磁南北方向的磁子午线上的一条极光强度数据，以时间为序列，排列成极光活动图.Keogram图能够大致反映极光图像亮度的变化，但是它只提取了每幅图像中的一条线的信息，反映的极光活动是有限的.变化曲线是度量极光活动与之前发生的极光活动的变化情况，曲线上的值越大表示极光的变化越剧烈，值越小说明极光的运动模式发生的变化越小.

实验结果显示，变化曲线完全能够准确定位极光事件的发生，与人工标记的结果非常吻合.图 4a是用于检测2003-12-24 11∶00—15∶00(世界时)极光活动的变化曲线，图 4b是对应时段的Keogram图.从图 4a的变化曲线可以看出，极光在这一时段发生了5次重大变化.第一个局部峰值位于12∶23(世界时)，图 4b的Keogram图显示：在2003-12-23 12∶23(世界时)以前，极光主要发生在黄河站以北的区域(全天空极光图像的上半部分)，而且极光的亮度较暗以及运动速度相对较低；在12∶23(世界时)以后，极光迅速开始向南移动.直到第二个局部峰值(12∶55 世界时)，这一运动趋势逐渐消退，极光的亮度开始变暗.在12∶55(世界时)之后，极光的活动发生了明显的变化，明亮的极光又一次出现，并且向极向方向移动.第三个局部峰值发生在13∶07(世界时)，对照Keogram图看出正是在这一时刻，之前的极光极向运动趋势结束，极光突然消失.从13∶23(世界时)开始，极光亮度突然变亮，而且极光具有极向运动趋势.这一趋势一直延续到最后一个局部峰值(13∶53(世界时)).在13∶53(世界时)之后，极光逐渐地消失.以上对比可见，这一时段极光的亮度、形状以及发生位置都发生了很大的变化，变化曲线准确地检测到极光活动变化的发生.

图 4

Fig. 4

图 4 2003-12-24 11∶00—15∶00 (世界时) 极光活动的变化曲线与对应时段的Keogram图Fig. 4 The change curve of auroral activities at 11∶00—15∶00 (UT) on 2003-12-24 and the corresponding Keogram during this interval

图 5是本文方法用于检测2003-12-25 12∶00—15∶00(世界时)极光活动的变化曲线与Keogram图的比较.图 5a的变化曲线显示这段极光活动发生了5次重大变化.从图 5b可以看出这一时段的Keogram图连贯性较低，说明极光活动的变化非常突然.对比图 5(a，b)可见，变化曲线的5个局部峰值发生的位置准确的定位了极光发生突变的时刻.不仅如此，变化曲线值的大小直接反映了极光活动变化的剧烈程度.例如，在 2003-12-25 12∶44(世界时)，变化曲线的值达到了最大值，而且变化曲线的波形也比较尖锐，说明极光活动在这个时刻发生了非常明显的突变.图 5b的Keogram图显示在12∶44(世界时)之前，天空中完全没有极光发生，在12∶44(世界时)之后，天空出现了非常明亮的极光现象，并且极光朝着极向方向运动.相比之下，出现在变化曲线上13∶51(世界时)的局部峰值相对较低，变化曲线的波形也比较缓和，说明极光在此时发生的变化并不十分剧烈.对照Keogram图也可以看出，虽然13∶51(世界时)时极光的亮度和所处的位置都发生了变化，但是变化相对比较缓和.

图 5

Fig. 5

图 5 2003-12-25 12∶00—15∶00 (世界时) 极光活动的变化曲线与对应时段的Keogram图的比较Fig. 5 The change curve of auroral activities at 12∶00—15∶00 (UT) on 2003-12-25 and the corresponding Keogram during this interval

图 6a是采用变化曲线检测到的发生在2004-01-16 04∶50—05∶30(世界时)的极光活动变化，图 6b是对应时段的Keogram图.由图 6a的变化曲线可以看出极光在这一时段的活动十分活跃，而图 6b的Keogram图未能非常明显的显示出这些极光活动的变化.按照变化曲线的指示，通过观察极光序列我们发现，发生在这一时段的极光形态、亮度、位置以及运动模式变化非常迅速.变化曲线出现最大峰值的位置在05∶07∶53(世界时)，图 6(c，d)显示的是发生在这个时刻前后的全天空极光图像.由所示的全天空极光图像序列可以看出，在05∶07∶53(世界时)之前，极光出现了明显的涡旋结构，在05∶07∶53(世界时)之后涡旋结构完全消失.观察图 6b的Keogram图在05∶07∶53时刻虽然也发生了变化，但是不足以引起人们的注意.所以，对于连续涡旋结构这种重要的极光现象，通过浏览Keogram图则无法检测到它的发生.

图 6

Fig. 6

图 6 2004-01-16 04∶50—05∶30(世界时)极光活动的变化曲线以及对应时段的Keogram图和ASI图像Fig. 6 The change curve of auroral activities at 04∶50—05∶30 (UT) on 2004-01-16, as well as the corresponding Keogram and ASI images during this interval

对没有出现在南北子午线上的亮度模式，通过观察Keogram图也无法发现这些极光活动.图 7a是采用本文方法检测到的发生在2004-01-17 05∶20—06∶00(世界时)的极光活动变化，图 7b的Keogram图显示的是同一时段子午线上极光亮度的变化情况.由变化曲线可以看出极光在05∶32(世界时)附近发生了突变，这一结果无法从Keogram图看出，但是可以由图 7(c，d)所示的全天空极光图像序列得到证实.在05∶32∶02(世界时)之前，天空中出现的主要亮班区域在图像的右上角，斑状和射线结构比较明显，经过四分钟的发展，这一极光亮斑在05∶32∶02(世界时)逐渐消失.从05∶32∶02(世界时)开始出现弧状结构，一条明亮的极光弧从左边开始出现，之后逐渐生长直至消失，之后又有一条具有褶皱结构的极光弧从天顶位置开始发展.由于Keogram图是抽取每一帧极光图像中沿地磁南北方向的磁子午线上的一条极光强度数据，对没有出现在南北子午线上的亮度模式，通过观察Keogram图则无法发现这一极光活动.综上所述，如果通过观察Keogram图来浏览极光变化，极有可能会错过许多重要细节和二维信息.

图 7

Fig. 7

图 7 2004-01-17 05∶20—06∶00(世界时)极光活动的变化曲线以及对应时段的Keogram图和ASI图像Fig. 7 The change curve of auroral activities at 05∶20—06∶00 (UT) on 2004-01-17, as well as the corresponding Keogram and ASI images during this interval

面对日积月累产生的海量极光数据，本文提出了一种基于全天空极光运动场的自动检测极光活动变化的方法.使用这一方法用于检测2003—2004年北极黄河站的越冬观测数据，实验结果显示，本文的方法能够有效的检测到极光活动的变化，准确定位了极光发生突变的位置.因为Keogram图只是提取了极光图像上一条线的信息，它无法表现极光丰富多彩的二维形态特征和运动模式.与传统的Keogram图相比，变化曲线反映的是图像序列形态、位置、速度等更丰富信息的变化情况.如果能够结合Keogram图和变化曲线将能够更快速、有效、全面地显示极光活动的变化.

经分析，我们认为本文提出的方法之所以能够有效检测极光活动，主要有以下两个原因：1)本方法考虑到了极光的流体性质，针对极光形态多变、自发光的特点，引入基于连续方程作为数据约束的方法，克服了现有方法采用的亮度不变假设与极光现象不符合的问题，提高了极光运动场估计的准确性.2)极光的运动具有不确定的时间性和空间性，但是邻域内像素的运动具有明显的相关性.本文基于提取出的极光运动场，采用基于局部向量差的时空统计表征极光序列，局部向量差考虑了邻域内的局部相对运动场.本文的研究成果不但实现了快速浏览极光活动和准确定位极光事件的发生，而且还为从统计学的角度分析极光活动周期奠定了基础.在今后的研究中，我们将就此问题展开研究.