CN105631471A - 融合单帧特征和动态纹理模型的极光序列分类方法 - Google Patents

Abstract

本发明公开了一种融合单帧特征和动态纹理模型的极光序列分类方法。主要解决现有技术特征提取不全面，对极光序列分类准确率不高的问题。实现步骤为：1.获取训练极光序列和测试极光序列；2.提取测试极光序列的单帧词袋BOW-SIFT特征，组成特征序列；3.根据特征序列分别提取测试极光序列的特征观测矩阵和特征状态转移矩阵；4.根据特征观测矩阵和特征状态转移矩阵求测试极光序列到训练极光序列的马丁距离；5.根据马丁距离，对测试极光序列进行最近距离分类。本发明能实现对四类极光序列的计算机自动分类，且具有分类准确率高的优点，可用于图像序列的特征提取和识别。

Description

融合单帧特征和动态纹理模型的极光序列分类方法

技术领域

本发明属于图像处理技术领域，特别涉及四种形态极光序列的机分类方法，可用于图像序列的特征提取和计算机图像识别。

背景技术

极光是太阳风通过日侧极隙区注入到地球磁层时，沉降粒子沿磁力线与地球高层大气相互作用而产生的绚丽光辉。极光是极区空间天气物理过程的观测窗口，直观地反映了太阳风与地磁层的耦合过程，蕴含着大量日地空间的电磁活动信息，有着深刻的研究意义。

中国北极黄河站的全天空成像系统(All-skyCamera)同时对极光的三个典型谱段427.8nm，557.7nm和630.0nm进行连续观测，产生数以万计的极光图像，数据量巨大。WangQ等人在文章“Spatialtexturebasedautomaticclassificationofdaysideaurorainall-skyimages.JournalofAtmosphericandSolar-TerrestrialPhysics,2010,72(5):498-508.”中将极光按形态分为弧状、辐射状、帷幔状和热点状四类，并得出了四种极光类型的统计分布规律。Pedersen等人发表的文章“PedersenTR,GerkenEA.Creationofvisibleartificialopticalemissionintheaurorabyhigh-powerradiowaves.Nature,2005,433(7025):498-500”，胡泽俊等人发表文章“HuZJ,YangH,HuangD,etal.Synopticdistributionofdaysideaurora:Multiple-wavelengthall-skyobservationatYellowRiverStationinNy-Alesund,Svalbard.J.Atmos.Sol.-Terr.Phys.,2009,71(89):794-804”以及文章“LorentzenDA,MoenJ,OksavikK,etal.Insitumeasurementofanewlycreatedpolarcappatch.J.Geophys.Res.,2010,115(A12).”，提供了大量的研究材料，证明不同形态的极光对应着不同的磁层边界层动力学过程。如何准确高效地对极光视频进行分类，既是揭示其磁层源区动力学过程的一个关键，也是其发生机制研究的重要环节，然而极光形态和动态变化复杂，无疑为极地科研工作者带来了巨大的困难.

计算机图像识别和分析技术的发展为海量极光数据分类研究提供了可能。2004年，等人在文章“SyrjasuoM,PartamiesN.Numericimagefeaturesfordetectionofaurora[J].GeoscienceandRemoteSensingLetters,IEEE,2012,9(2):176-179.”中将计算机视觉的方法引入极光图像自动分类系统中，该方法从分割后的极光区域中提取Fourier算子作为特征，通过最近邻方法实现极光图像的自动分类，由于受到分割算法的影响，该方法仅对形状特征明显的弧状极光图像分类效果良好；Wang等人于2007年在文章“WangQian,LiangJimin,HuZeJun,HuHaiHong,ZhaoHeng,HuHongQiao,GaoXinbo,YangHuigen.Spatialtexturebasedautomaticclassificationofdaysideaurorainall-skyimages.JournalofAtmosphericandSolar-TerrestrialPhysics,2010,72(5):498–508.”中使用主成分分析法PCA对极光图像的灰度特征进行提取，提出了一种基于表象的极光分类方法，在冕状极光分类研究方向取得了一定进展；2008年，Gao等人发表文章“L.Gao,X.B.Gao,andJ.M.Liang.DaysidecoronaautoradetectionbasedonsampleselectionandadaBoostalgorithm.J.ImageGraph,2010,15(1):116-121.”，提出基于Gabor变换的极光图像分类方法，采用了局部Gabor滤波器提取图像特征，在确保计算精度的情况下降低了特征冗余信息，取得了较好的分类效果；2009年，Fu等人在文章“FuRu,JieLiandX.B.Gao..Automaticauroraimagesclassificationalgorithmbasedonseparatedtexture.Proc.Int.Conf.RoboticsandBiomimetics,2009:1331-1335.”中将形态学成分分析(MCA)与极光图像处理相结合，从经过MCA分离后所得到的极光纹理子图中提取特征，用于弧冕两类极光图像的分类，提高了弧冕极光分类的正确率。后续的相关研究还有：Han等人在文章“BingHan,XiaojingZhao,DachengTao,etal.DaysideauroraclassificationviaBIFs-basedsparserepresentationusingmanifoldlearning.InternationalJournalofComputerMathematics.Publishedonline:12Nov2013.”中又提出基于BIFs特征和C均值分类的极光分类法；Yang等人在文章“YangXi,LiJie,HanBing,GaoXinbo.Wavelethierarchicalmodelforauroraimagesclassification.JournalofXidianUniversity,2013,40(2):18-24.”中提出多层小波变换来表示极光图像特征，取得了较高的分类准确率；2013年，Han等人在文章“HanB,YangC,GaoXB.AuroraimageclassificationbasedonLDAcombiningwithsaliencyinformation.RuanJianXueBao/JournalofSoftware,2013,24(11):2758-2766.”中引入隐含狄利克雷分布模型LDA，并结合图像显著性信息，又进一步提高了极光图像的分类准确率.

然而，关于极光序列的分类的相关工作仍然还比较少。相关的进展主要有：Yang等人于2013年在文章“YangQiuju.AuroralEventsDetectionandAnalysisBasedonASIandUVIImages[D].Xi’an:Xidianuniversity,2013.”中提出基于隐马尔科夫模型的极光序列分类理论，但是该方法本质仍是基于单幅图像特征；另外，Han在文章“HanB,LiaoQ,GaoXB.Spatial-Temporalpolewardvolumelocalbinarypatternsforaurorasequenceseventdetection.RuanJianXueBao/JournalofSoftware,2014,25(9):2172-2179.”中构造了空时极向LBP表征算子提取极光序列的STP-LBP特征，用于检测极光视频中的极向运动现象。但是，目前这些对于极光序列分类主要是基于单帧特征和动态特征中的一种，特征提取方法单一，导致分类准确率和分类效率都比较低。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种融合单帧特征和动态纹理模型的极光序列分类方法，以综合考虑极光序列的单帧特征和动态变化特征，对极光的四类形态序列进行完整的特征提取，提高极光序列的分类准确率。

实现本发明上述目的的技术方案是：提取极光序列的单帧BOW-SIFT特征，组成特征序列；通过SVD分解提取特征序列的动态纹理模型参数；然后根据动态纹理模型参数构建马丁距离，用最近距邻分类器完成对极光序列的分类，其实现步骤包括如下：

(1)从已标记类别的极光序列数据库中任取N个序列组成训练集{S₁,S₂,...,S_j,...,S_N}，S_j是第j个训练集样本，j＝1,2,...,N，将剩余极光序列组成测试集，从测试集中取一个样本作为测试极光序列St，其中St＝{I₁,I₂,...,I_t,...,I_τ}，I_t表示极光序列St的第t帧图像，t＝1,2,...,τ，τ表示序列长度；

(2)提取极光序列St中任意帧图像I_t的单帧BOW-SIFT特征y_t，组成特征序列F＝{y₁,y₂,...,y_t,...,y_τ}；

(3)根据特征序列F提取极光序列St的特征观测矩阵C和特征动态纹理矩阵X；

(4)用得到的特征动态纹理矩阵X，计算极光序列St的特征状态转移矩阵A，

$\begin{array}{c}A=\arg \min \left|\right|X_{2,\mathrm{...},\mathrm{\τ}}-{\mathrm{AX}}_{1,\mathrm{...},\mathrm{\τ}-1}|{|}_F^2\\ =\left(X_{2,\mathrm{...},\mathrm{\τ}}{X}_{1,\mathrm{...},\mathrm{\τ}-1}\right){\left(X_{1,\mathrm{...},\mathrm{\τ}-1}{X}_{1,\mathrm{...},\mathrm{\τ}-1}^T\right)}^{-1}\end{array},$

式中X_1,...,τ-1表示特征动态纹理矩阵X的第1列到第τ-1列所组成的矩阵；X_2,...,τ表示矩阵X的第2列到第τ列所组成的矩阵；表示求矩阵的F范数；

(5)根据特征观测矩阵和特征状态转移矩阵，求测试极光序列St到训练集样本S_j的马丁距离d²(St,S_j)；

(6)将步骤(5)中得到的N个马丁距离d²(St,S_j)按从小到大的顺序排列，取出最小的马丁距离对应的极光序列S_min，并将与该极光序列S_min形态相似的测试极光序列St分为与极光序列S_min同一类，完成对测试极光序列St的分类。

本发明具有如下优点：

本发明由于采用先提取极光序列单帧特征，组成特征序列；然后利用SVD分解提取极光特征序列的动态纹理特征，充分融合了极光序列的单帧特征和动态变化特征，克服了现有技术只分析单帧特征和动态特征中的一种，造成特征提取不全面，极光序列分类准确率不高的缺点，使得该方法能够提高对四类形态极光序列的准确率。

附图说明

图1是本发明的实施例原理图。

具体实施方式

以下结合附图对本发明的内容和效果进行进一步描述。

参照图1，本发明的实现方案如下：

步骤1：获取训练极光序列和测试极光序列。

从已标记类别的极光序列数据库中任取N个序列组成训练集{S₁,S₂,...,S_j,...,S_N}，S_j是第j个训练集样本，j＝1,2,...,N，将剩余极光序列组成测试集；

从测试集中取一个样本作为测试极光序列St，其中St＝{I₁,I₂,...,I_t,...,I_τ}，I_t表示极光序列St的第t帧图像，t＝1,2,...,τ，τ表示序列长度。

步骤2：提取极光序列St中任意帧图像I_t的单帧词袋BOW-SIFT特征y_t，组成特征序列F。

(2a)将512×512大小的图像I_t进行滑窗网格划分，滑窗大小32×32，每次滑动16像素，得到961个滑窗区域；

(2b)将(2a)得到的全部961个32×32区域均提取128维的尺度不变特征变换SIFT特征，得到图像I_t的尺度不变特征变换SIFT特征矩阵F_S，F_S∈R^916×128；

(2c)将训练集内所有帧图像的尺度不变特征变换SIFT特征矩阵F_S进行K-means聚类，得到K个聚类中心Fc＝{Fc₁,Fc₂,...,Fc_d,...,Fc_K}，将Fc称为字典，将Fc_d称为第d个视觉单词，其中d＝1,2,...,K，K为字典大小；

(2d)用字典Fc对尺度不变特征变换SIFT特征矩阵F_S进行量化表示，得到图像I_t的词袋BOW-SIFT特征y_t：

(2d1)将SIFT特征矩阵F_S进行转置，并表示为F_S＝{fs₁,fs₂,...,fs_z,...,fs_W}，其中，fs_z∈R¹²⁸，z＝1,2,...,W，W＝916；

(2d2)分别计算fs_z到视觉单词Fc_d的欧式距离，并用与fs_z距离最小的视觉单词Fc_d替换特征矩阵F_S中的fs_z，得到替换后的特征矩阵F_R；

(2d3)统计替换后的特征矩阵F_R中每个视觉单词Fc_d出现的频率，得到长度为K的词频向量y_t，称为图像I_t的词袋BOW-SIFT特征。

步骤3：提取极光特征序列F的特征观测矩阵C和特征动态纹理矩阵。

根据动态纹理特征模型理论，求序列的特征观测矩阵和特征动态纹理矩阵技术包括：子空间参数鉴别N4SID法、输出状态误差MOESP法，以及基于最大似然估计原理的SVD闭式求解法等。为提高运算效率，本实例采用SVD闭式求解法，其具体操作步骤描述如下：

(3a)用特征序列F减去所有特征帧的均值得到归一化特征矩阵Y：

$Y=\{y\left(1\right)-\stackrel{\‾}{y},y\left(2\right)-\stackrel{\‾}{y},...,y\left(t\right)-\stackrel{\‾}{y},...,y\left(\mathrm{\τ}\right)-\stackrel{\‾}{y}\};$

(3b)对归一化特征矩阵Y进行SVD分解，即Y＝USV^T，其中，U为左基矩阵，S为奇异值矩阵，V为右基矩阵，T为求矩阵转置的运算符号；

(3c)用左基矩阵U表示极光特征序列的特征观测矩阵C，即C＝U；用奇异值矩阵S和右基矩阵V求极光特征序列的特征动态纹理矩阵X＝SV^T，其中X∈R^K×τ。

步骤4：用得到的特征动态纹理矩阵X，计算特征状态转移矩阵A：

$\begin{array}{c}A=\arg \min \left|\right|X_{2,\mathrm{...},\mathrm{\τ}}-{\mathrm{AX}}_{1,\mathrm{...},\mathrm{\τ}-1}|{|}_F^2\\ =\left(X_{2,\mathrm{...},\mathrm{\τ}}{X}_{1,\mathrm{...},\mathrm{\τ}-1}\right){\left(X_{1,\mathrm{...},\mathrm{\τ}-1}{X}_{1,\mathrm{...},\mathrm{\τ}-1}^T\right)}^{-1}\end{array},$

式中X_1,...,τ-1表示特征动态纹理矩阵X的第1列到第τ-1列所组成的矩阵；X_2,...,τ表示矩阵X的第2列到第τ列所组成的矩阵；表示求矩阵的F范数。

步骤5：根据特征观测矩阵和特征状态转移矩阵，求测试极光序列St到训练集样本S_j的马丁距离d²(St,S_j)。

(5a)计算训练集{S₁,S₂,...,S_j,...,S_N}中训练集样本S_j的训练特征观测矩阵C_j和训练特征状态转移矩阵A_j，得到训练特征观测矩阵集{A₁,A₂,...,A_j,...,A_N}和训练特征状态转移矩阵集{C₁,C₂,...,C_j,...,C_N},其中j＝1,2,...,N；

(5b)用测试极光序列St的特征观测矩阵C和特征状态转移矩阵A以及训练集样本的训练特征观测矩阵C_j和训练特征状态转移矩阵A_j构建扩展观测矩阵：

$\begin{array}{c}\mathrm{\Φ}t=\{C^T,A^T{C}^T,{\left(A^T\right)}^2{C}^T,\mathrm{...},{\left(A^T\right)}^i{C}^T,\mathrm{...},{\left(A^T\right)}^B{C}^T\}\\ {\mathrm{\Φ}}_j=\{C_j^T,A_j^T{C}_j^T,{\left(A_j^T\right)}^2{C}_j^T,\mathrm{...},{\left(A_j^T\right)}^i{C}_j^T,\mathrm{...},{\left(A_j^T\right)}^B{C}_j^T\}\end{array};$

其中，Φt为测试扩展观测矩阵，Φ_j为训练扩展观测矩阵，T为求矩阵转置的运算符号，i＝1,2,...,B，B表示截断常数，其取值为一个大于5的常数。

(5c)计算测试扩展观测矩阵Φt到训练扩展观测矩阵Φ_j的第i个特征角θ_i(Φt,Φ_j)，即：

${\mathrm{\θ}}_i(\mathrm{\Φ}t,{\mathrm{\Φ}}_j)=\arccos \left(\frac{\left|{\left(A^T\right)}^i{C}^T{\left(A_j^T\right)}^i{C}_j^T\right|}{\left|{\left(A^T\right)}^i{C}^T\left|\right|{\left(A_j^T\right)}^i{C}_j^T\right|}\right),$

其中，|·|表示求矩阵的模值，arccos表示反余弦运算；

(5d)根据特征角θ_i(Φt,Φ_j)计算测试极光序列St到训练集样本S_j的马丁距离:

$d^2(St,S_j)=-log\underset{i=1}{\overset{B}{\Π}}{\cos }^2{\mathrm{\θ}}_i(\mathrm{\Φ}t,{\mathrm{\Φ}}_j).$

步骤6：将步骤(5)中得到的N个马丁距离d²(St,S_j)按从小到大的顺序排列，取出最小的马丁距离对应的极光序列S_min，并将与该极光序列S_min形态相似的测试极光序列St分为与极光序列S_min同一类，完成对测试极光序列St的分类。

本发明的效果通过以下仿真实验进一步说明：

1.仿真条件与方法：

硬件平台为：IntelCorei5、2.93GHz、3.45GBRAM；

软件平台为：Windows7操作系统下的MATLABR2012b；

实验的数据：对我国北极黄河站的全天空极光数据中2003年12月到2004年共115557幅G波段图像进行手工标记，去除天气等因素造成的无效数据，标记出93个辐射状极光序列，102个弧状极光序列，73个热点状极光序列，95个帷幔状极光序列，序列长度为15帧到35帧之间，组成典型极光序列数据库进行四分类实验。

2.仿真内容与结果：当训练集序列数N＝30,50,80,100,120个时，分别使用现有的动态纹理模型和本发明所提出方法对极光序列进行分类，得到如下表所示的分类准确率对比。

表1极光序列分类准确率

从表1可以看出，采用本发明所提出的融合单帧特征和动态纹理模型的极光序列分类方法有效提高了极光序列的分类准确率。

Claims (4)

1.一种融合单帧特征和动态纹理模型的极光序列分类方法，包括如下步骤：

(1)从已标记类别的极光序列数据库中任取N个序列组成训练集{S₁,S₂,...,S_j,...,S_N}，S_j是第j个训练集样本，j＝1,2,...,N，将剩余极光序列组成测试集，从测试集中取一个样本作为测试极光序列St，其中St＝{I₁,I₂,...,I_t,...,I_τ}，I_t表示极光序列St的第t帧图像，t＝1,2,...,τ，τ表示序列长度；

(2)提取极光序列St中任意帧图像I_t的单帧词袋BOW-SIFT特征y_t，组成特征序列F＝{y₁,y₂,...,y_t,...,y_τ}；

(3)根据特征序列F提取极光序列St的特征观测矩阵C和特征动态纹理矩阵X；

(4)用得到的特征动态纹理矩阵X，计算极光序列St的特征状态转移矩阵A，

$\begin{array}{c}A=\arg \min \left|\right|X_{2,...,\mathrm{\τ}}-{\mathrm{AX}}_{1,...,\mathrm{\τ}-1}|{|}_F^2\\ =\left(X_{2,...,\mathrm{\τ}}{X}_{1,...,\mathrm{\τ}-1}\right){\left(X_{1,...,\mathrm{\τ}-1}{X}_{1,...,\mathrm{\τ}-1}^T\right)}^{-1}\end{array},$

式中X_1,...,τ-1表示特征动态纹理矩阵X的第1列到第τ-1列所组成的矩阵；X_2,...,τ表示矩阵X的第2列到第τ列所组成的矩阵；表示求矩阵的F范数；

(5)根据特征观测矩阵和特征状态转移矩阵，求测试极光序列St到训练集样本S_j的马丁距离d²(St,S_j)；

(6)将步骤(5)中得到的N个马丁距离d²(St,S_j)按从小到大的顺序排列，取出最小的马丁距离对应的极光序列S_min，并将与该极光序列S_min形态相似的测试极光序列St分为与极光序列S_min同一类，完成对测试极光序列St的分类。

2.根据权利要求1所述的方法，其中所述步骤(2)中提取图像I_t的单帧BOW-SIFT特征y_t，按如下步骤进行：

(2a)将512×512大小的图像I_t进行滑窗网格划分，滑窗大小32×32，每次滑动16像素，得到961个滑窗区域；

(2b)将(2a)得到的全部961个32×32区域均提取128维的SIFT特征，得到图像I_t的尺度不变特征变换SIFT特征矩阵F_S，F_S∈R^916×128；

(2c)将训练集内所有帧图像的尺度不变特征变换SIFT特征矩阵F_S进行K-means聚类，得到K个聚类中心Fc＝{Fc₁,Fc₂,...,Fc_d,...,Fc_K}，将Fc称为字典，将Fc_d称为第d个视觉单词，其中d＝1,2,...,K，K为字典大小；

(2d)用字典Fc对尺度不变特征变换SIFT特征矩阵F_S进行量化表示，得到图像I_t的词袋BOW-SIFT特征y_t：

(2d1)将SIFT特征矩阵F_S进行转置，并表示为F_S＝{fs₁,fs₂,...,fs_z,...,fs_W}，其中，fs_z∈R¹²⁸，z＝1,2,...,W，W＝916；

(2d2)分别计算fs_z到视觉单词Fc_d的欧式距离，并用与fs_z距离最小的视觉单词Fc_d替换特征矩阵F_S中的fs_z，得到替换后的特征矩阵F_R；

(2d3)统计替换后的特征矩阵F_R中每个视觉单词Fc_d出现的频率，得到长度为K的词频向量y_t，称为图像I_t的词袋BOW-SIFT特征。

3.根据权利要求1所述的方法，其中所述步骤(3)根据特征序列F提取极光序列St的特征观测矩阵C和特征动态纹理矩阵X，按如下步骤进行：

(3a)用特征序列F减去所有特征帧的均值得到归一化特征矩阵Y：

$Y=\{y\left(1\right)-\stackrel{\‾}{y},y\left(2\right)-\stackrel{\‾}{y},...,y\left(t\right)-\stackrel{\‾}{y},...,y\left(\mathrm{\τ}\right)-\stackrel{\‾}{y}\};$

(3b)对归一化特征矩阵Y进行SVD分解，即Y＝USV^T，其中，U为左基矩阵，S为奇异值矩阵，V为右基矩阵，T为求矩阵转置的运算符号；

(3c)用左基矩阵U表示极光特征序列的特征观测矩阵C，即C＝U；用奇异值矩阵S和右基矩阵V求极光特征序列的特征动态纹理矩阵X＝SV^T，其中X∈R^K×τ。

4.根据权利要求1所述的方法，其中所述步骤(5)中根据特征观测矩阵和特征状态转移矩阵，求解测试极光序列St到训练集样本S_j的马丁距离d²(St,S_j)，按如下步骤进行：

(5a)计算训练集{S₁,S₂,...,S_j,...,S_N}中训练集样本S_j的训练特征观测矩阵C_j和训练特征状态转移矩阵A_j，得到训练特征观测矩阵集{A₁,A₂,...,A_j,...,A_N}和训练特征状态转移矩阵集{C₁,C₂,...,C_j,...,C_N},其中j＝1,2,...,N；

(5b)用测试极光序列St的特征观测矩阵C和特征状态转移矩阵A以及训练集样本的训练特征观测矩阵C_j和训练特征状态转移矩阵A_j构建扩展观测矩阵：

Φt＝{C^T,A^TC^T,(A^T)²C^T,...,(A^T)ⁱC^T,...,(A^T)^BC^T}

${\mathrm{\Φ}}_j=\{C_j^T,A_j^T{C}_j^T,{\left(A_j^T\right)}^2{C}_j^T,...,{\left(A_j^T\right)}^i{C}_j^T,...,{\left(A_j^T\right)}^B{C}_j^T\};$

其中，Φt为测试扩展观测矩阵，Φ_j为训练扩展观测矩阵，T为求矩阵转置的运算符号，i＝1,2,...,B，B表示截断常数，其取值为一个大于5的常数。

(5c)提取测试扩展观测矩阵Φt到训练扩展观测矩阵Φ_j的第i个特征角θ_i(Φt,Φ_j)，即：

${\mathrm{\θ}}_i\left(\mathrm{\Φ}t,{\mathrm{\Φ}}_j\right)=\arccos \left(\frac{\left|{\left(A^T\right)}^i{C}^T{\left(A_j^T\right)}^i{C}_j^T\right|}{\left|{\left(A^T\right)}^i{C}^T\right|\left|{\left(A_j^T\right)}^i{C}_j^T\right|}\right),$

其中，|·|表示求矩阵的模值，arccos表示反余弦运算；

(5d)根据特征角θ_i(Φt,Φ_j)计算测试极光序列St到训练集样本S_j的马丁距离:

$d^2(St,S_j)=-log\underset{i=1}{\overset{B}{\Π}}{\cos }^2{\mathrm{\θ}}_i(\mathrm{\Φ}t,{\mathrm{\Φ}}_j).$