基于支持向量机的云图自动识别和提取方法∗

王利文 贾 鹏† 蔡冬梅 刘慧根

(1太原理工大学物理与光电工程学院太原030024)

(2南京大学天文与空间科学学院南京210034)

1 引言

时域天文学是目前非常活跃的一个研究领域,其研究对象主要包括超新星、变星及系外行星等光度连续变化的目标源和太阳系内的近地天体等快速运动目标.虽然这些天体位置与光度变化的时间尺度各不相同,且对于观测提出的具体要求也各有不同,但是总体上都要求在各自的时间窗口内保证观测数据的连续获取,因此大多使用多台地理上分散分布的中小口径望远镜,通过自动控制系统对天空进行连续成像观测[1−7].但这种观测手段给时域天文的观测和数据处理带来了许多新的挑战和问题[8]:在时域天文观测时,除了近地快速移动目标为线源外,大部分感兴趣的目标属于点源,在仪器端获取的大量的天文数据中,会有部分数据受到天空中云层的污染.在对这些数据中的天文目标进行探测和光度测量时,云层往往会影响观测精度,甚至会完全干扰观测[9].如图1所示为仪器端获取的图像数据,其中(a)–(c)为没有云的正常图像,(d)–(f)为受云干扰的图像.

图1 图像示例.(a)–(c)正常图像;(d)–(f)含云图像Fig.1 Image examples.(a)–(c)normal images;(d)–(f)images with cloud

由于云的影响,光透过云层强度会发生衰减,给测光带来误差;同时来自地面的光会从云区直接反射到成像系统,使有云区比无云区更亮一些,严重影响对暗弱暂现源的提取.基于此,为了给测光和暗弱暂现源提取时提供参考,有必要首先对图像中的云进行提取并根据云的灰度变化建立指标图.指标图是一幅与原始数据大小相等的图片,其每个像元的灰度信息可以反映出获取的数据中的云的轮廓和云的灰度信息.但是如果图像中没有云,获得的指标图对于观测将变得毫无意义.因此,需要在提取云之前筛选出含云的图像.由于时域天文观测的数据量大,只依靠人工筛选耗时且繁重,所以有必要构建云图快速分类系统.之后根据分类系统进一步提取图像中的云层,并建立指标图以便于科学研究[10−11].

但是,图像的数据特性评价一直是一个非常困难的问题,对于天文图像尤甚.由于外界干扰因素很多,且不同的干扰因素对于不同天文观测任务的影响也各不相同.因此,从成像的物理过程分析,通过确定几个经验参数及其分布范围来直接判断图像中是否存在云非常困难.除了从物理过程直接分析外,另一种思路就是从数据特征出发,将存在云的天文图像进行标记,结合云图像的形态特征,通过机器学习的方法建立符合数据特征的分类器.

在星系分类和光谱识别中[12−13],支持向量机(Support Vector Machine,SVM)是一种广泛应用的分类器.该分类器的基本策略是保证不同类别的数据具有最大的分类间隔.由于这类求间隔最大化的问题往往可以转化为凸二次规划问题,因此与神经网络、随机森林和决策树等工具相比,SVM可以在数据量较少的情况下快速得到需要的分类器[14−15].这一特性降低了数据积累的要求,同时减少了人工设置标签的工作量.此外,当数据在低维空间不可分时,通过核函数映射可以将分类问题转化为求取高维空间分离超平面的问题.本文将在第2章尝试从云图像的特征出发,通过分析云图像的特征,确定特征描述量构建分类系统,并进行了实验测试;第3章对分类器获得的有云的图像进行了指标图提取;最后对文章进行了总结.

2 天文云图快速分类系统

由于云是一种具有2维结构的图像,我们需要从它的特征出发,利用一些参数对其进行描述.一般来说,云层具有一定的延展特征和纹理特性,且云层的厚薄会导致灰度差异,为此,我们把图像的灰度不一致度和纹理特征设定为天文图像云层的评价指标,采用灰度不一致度和纹理特征共5个指标共同对云进行识别.

2.1 云图特征提取

2.1.1 灰度不一致度

图像的灰度不一致度是反映图像背景灰度的一个重要指标,作为天文图像的一个评价指标,图像的灰度不一致度G的表达式如下:

式中,Ps和Pn分别为用局部图像均方差的最大值与最小值.根据实际观测时不同设备的工作情况可以得到典型的星像尺度,再根据其大小调节模板尺寸,本文选用9×9的模板遍历全图.

根据上述图像的灰度不一致度的定义,我们可以发现该特征量所反映的是图像在特征尺度内,灰度值整体波动(背景)和局部特殊值(恒星或其他天文目标)的不一致度.图像的灰度不一致度的高低对图像质量有着重要的参考价值,当有云存在时,图像的灰度不一致度会发生剧烈变化.

2.1.2 纹理特征

在图像灰度不一致度的基础上,采用统计信息——纹理特征可以更好地对云进一步描述.灰度共生矩阵能够反映图像灰度关于方向、相邻间隔、变化幅度的综合信息,所以其适用于对云图像进行纹理识别,进而可以用于检测天文图像中是否有云的存在.满足一定空间关系的灰度共生矩阵元素为:

式中,i和j分别为矩阵g的行数和列数,M和N为矩阵g的行列数,g为图像f的灰度共生矩阵,#(x)表示集合中x元素的个数.若(x1,z1)与(x2,z2)为矩阵f中的两个点,两点之间距离为d,两者与坐标横轴的夹角为q,则可以得到各种间距及角度的灰度共生矩阵g(i,j,d,q).我们这里研究的对象是天文图像,主要的观测目标图像一般都具有旋转对称的特点,这一特点导致在不同方向的灰度共生矩阵对图像的纹理特征影响差异很小,所以这里选用d=1,q=0.考虑到望远镜自动观测时,获取的天文图像中云的特征,我们选用如下4个纹理特征对图像进行评价:

(1)能量能量即每个元素的平方和,是对图像平整度的衡量.如果g中的所有值均匀,则ASM的值较小;相反,在有云的情况下,部分区域有较大的值而其他区域值较小,则ASM的值较大,所以ASM可以用作对云图的一个识别指标.

其中,Mg和Ng分别为矩阵g的行列数.

(2)对比度对比度能够反映图像局部灰度变化的情况,云的存在会导致部分区域灰度变化很大,故有云图像的CON值较大,所以CON可以作为云图的另一个识别指标.

(3)逆差距逆差距反映了图像纹理的同质性,可以用来度量图像纹理局部变化的多少.由于云往往具有一定的延展结构,其纹理特性变化会比较大,所以IDM可以作为检测云图的又一重要指标.

(4)熵熵是图像所具有的信息量的度量,是对图像随机性的度量.当相似观测条件下的天文图像存在云时,熵值较不存在云时高.由于时域天文观测图像主要是点状目标,其熵值小于同样条件下包含云的图像.因此ENT可以作为检测云的最后一个重要指标.

2.2 云图分类器

2.2.1 SVM

SVM是由Vapnik等人提出的一种常用的分类算法[16].最早被用于解决如下的二分类问题:

其中,w=(w1,w2,···,wn)为决定超平面方向的法向量,n为x的维度,x表示数据点,b为超平面与原点之间的距离,SVM通过在高维数据集中寻找一个超平面来实现不同类别之间的数据几何间隔距离最大化.几何间隔s定义为:

其中,y表示这些数据点代表的类别,分别用1和–1表示,l为函数间隔,可以由下式求出:

则求几何间隔距离最大化问题可以转化为下面的凸优化问题:

上述问题可以使用二次规划优化包直接求解[17].但是,从上一小节的图像特征描述模型讨论中,我们发现各个评价指标之间并不是独立的,由于图像的复杂性,这些指标往往两两之间存在着关联.当云层较薄时,随着云层面积增加,灰度整体波动较大,ASM会加大,同时由于部分区域灰度变大,CON也会变化;但是,如果云层同时厚度增加,部分区域灰度不均匀性将降低,此时CON对应值将会不变,甚至变小.由于这些量关联性比较复杂,在可以满足需求的情况下,如果再增加其他参量,关联性可能会破坏特征空间的结构,导致SVM无法找到超平面,造成图像分类能力不佳.对于这一问题,SVM可以将线性特征映射到更高维特征空间来实现图像特征的可分[15].一个比较简单的办法就是通过核函数的技巧,将上述分类函数转化为:

其中,K(xi,xj)为核函数,ai为拉格朗日乘子,经过转化,凸优化问题成为了如下的形式:

2.2.2 SVM云图分类器构建过程

我们首先使用线性核作为SVM核函数,根据其对测试集的分类性能进行评价.当系统分类性能较差时(准确率小于90%),再进一步考虑将SVM核转换为非线性核.分类系统流程如图2所示.

从图2可以看出,整个云图像实时筛选系统包括如下步骤:

(1)读取每一幅原始图像的数据,并人工根据原始数据是否有云添加标签;

(2)将读取的数据分为两部分,一部分用作训练集,其余部分用作测试集.根据经验,训练集数目应该远大于特征维度,一般用100幅左右图像可以达到训练的目的;

(3)分别对训练集和测试集的5个维度进行特征提取,并对提取的特征归一化处理.特征空间维度如表1所示,其中Original feature是原始的特征值,Normalized feature为规范化的特征值,我们选X1作为参考,将所有特征值映射到和X1同一个量级,其中X2归一化为10000X1,X3归一化为Ln(X3),X4归一化为100X4,X5保持原始值不变;

(4)用训练集提取的特征量训练SVM分类器;

(5)用训练好的SVM分类器对测试集进行分类.X1–X5的归一化是为了消除不同特征量级差距的影响,使SVM训练时,保证数据在同一量级.需要注意的是:归一化是因数据而异的,也就是对于其他观测需要根据数据特性设计其他归一化方式,保证它们在量级上不会有差异[18].对图1中的示例图像,从左到右分别计算图像的特征值,计算的特征值如表2所示.

图2 系统流程图Fig.2 System flowchart

表1 特征空间的维度Table 1 Dimensions of feature space

2.3 GPU加速

实际工作时,该图像分类器将加载于控制计算机上对CCD采集的数据进行分类.但是在G提取中采用了全图遍历,使得运行的速度比较慢.处理大小为4096×4096的图像,对G这一特征提取就需要耗时300秒·幅−1,因此我们需要对算法进一步加速.

GPU又称图像处理器,具有强大的浮点运算和并行计算能力,与CPU相比更适合处理大量的并行数据.本文提出的办法在计算灰度不一致度时由于任意两次运算之间两两独立,所以适合采用GPU进行加速.望远镜拍摄的图片大小为M×N,遍历采用的模板大小为m×n,则对单幅图片需要运算的次数为(M−m+1)×(N−n+1).

经过测试,我们发现使用GPU加速后,对G这一特征的提取速度达到了0.43秒·幅−1,加速比接近700倍.

表2 正常图像与含云图像的特征值对比Table 2 The contrast of eigenvalues between normal images and images with cloud

2.4 分类实验结果与分析

2.4.1 实验数据

本文实验数据选用的是南京大学时域天文台(Time Domain Observatory,TIDO)18 cm望远镜在紫金山天文台盱眙观测站试观测时的数据.望远镜视场为12◦,使用安道尔IkonXL系列CCD,像元数为4096×4096.我们从中选取665幅图像作实验样本,对这些数据进行人工标记,其中有云的图片150幅,没有云的515幅.对图像有云、没云分别标注为类1、类2.

2.4.2 实验过程与结果

我们分别从类1、类2随机选取50%作为训练集,剩余的作为测试集数据.然后按照图2中所示的流程开始进行实验测试,最终获得的结果用图3的混淆矩阵表示出来.

图3 混淆矩阵Fig.3 The confusion matrix

从图3中可以看出:在332幅测试图像中,类1共有75幅,在实际测试中有1张图像被错分为类2,其错误率为1.3%;类2共有257幅,在实际测试中有3张被误分为类1,其错误率为1.2%.类1、类2的总体错误率为1.2%.这一分类准确率基本达到了我们的要求,为下一步提取云提供了参考.

3 云的提取

为了保证观测的空间和时间连续性,虽然含云的图像质量下降,但通过参考指标图,去除掉观测数据中含云的部分,剩余的数据仍然可以利用.图4为我们设计的分类算法识别出的有云的图像,下面以图4为例介绍云的指标图的提取过程.

图4 有云图像Fig.4 Image with clouds

由于云是整体具有一定轮廓和纹理的图像,从图像信息的频率域分析,其信息主要分布于中低频;从图像像元的灰度分布分析,当云的灰度值比较高时(比较厚),包含云的像元会在直方图上产生一个峰,当云的灰度值比较低时(比较薄),包含云的像元会扩展背景像元形成的峰的宽度.为此,我们设计了如下过程对图像中含云的部分进行提取:

(1)提取图像的灰度直方图,对图像进行高斯多峰拟合,并求拟合函数的拐点.由于图像的灰度值主要分布在图像的背景区域,图像灰度直方图近似高斯分布,故采用如下多高斯函数进行拟合:

其中,a1、b1、c1、a2、b2、c2、a3、b3、c3为拟合函数的系数, 求拟合函数系数,对h求2阶导数,与x轴的交点即为函数的拐点,分别取两个拐点值为背景灰度值的上下阈值,图5为灰度直方图的多高斯拟合曲线h.

图5 灰度直方图的多高斯拟合Fig.5 Multi-peaks Gaussian fitting of grayscale histogram

(2)分别取两个拐点值为背景灰度值的上下阈值,将这两个阈值之间的部分6等分,依次用不同的阈值绘制图像等高线并删除面积小于500的区域,用该阈值替代这一梯度的灰度值,所得的灰度梯度图如图6所示.

(3)对原始图像进行滤波,由于图像内背景和云是低频信息,而图像中的星像为高频信息,通过使用低通滤波器,可以降低个别亮星的干扰.为此我们选择高斯滤波,通过SExtractor提取星的位置信息,预估图像中星像的最大直径,结合经验选择合适大小的核对图像进行滤波,我们这里核大小选择为50.

(4)对滤波后的图像,按照(2)中所求的梯度,分别在该梯度进行如下操作:

所求的云的指标图如图7所示.从图中可以看到,大部分云的轮廓可以提取出来,能够为测光和后续暗弱目标提取的参考.但是亮星的影响比较难去除,是影响这个办法性能的主要因素.图8是更多的含云图像处理结果,其中(a)–(c)为含云图像,(d)–(f)为云的指标图.

图6 灰度梯度图Fig.6 The grayscale gradient of the image

图7 云的指标图Fig.7 The index figure of clouds

图8 含云图像及其处理结果.(a)–(c)含云图像;(d)–(f)云的指标图Fig.8 Images with colud and their process result.(a)–(c)images with colud;(d)–(f)the index figure of clouds

4 总结

本文根据时域天文观测数据处理需求,针对天文图像中云的特点提出了融合灰度不一致度与灰度纹理的图像特征,并据此构建了基于SVM的实时云图筛选系统,在筛选获得的图像基础上对图片中云的轮廓进行了提取.该系统的分类准确度可以达到98%以上,准确率高,鲁棒性强,同时也弥补了人工筛选中自动化程度低的不足.从实际应用角度,系统中采用了GPU技术,满足高速处理的需求.经过处理后获得的云的指标图可以大致反映云的轮廓,但是部分结果会受到亮星的影响,这一问题有待于进一步研究.

致谢感谢审稿人提出的建议和紫金山天文台孙荣煜对本研究的大力支持.