王利文 賈 鵬? 蔡冬梅 劉慧根

(1太原理工大學(xué)物理與光電工程學(xué)院太原030024)

(2南京大學(xué)天文與空間科學(xué)學(xué)院南京210034)

1 引言

時(shí)域天文學(xué)是目前非?；钴S的一個(gè)研究領(lǐng)域,其研究對(duì)象主要包括超新星、變星及系外行星等光度連續(xù)變化的目標(biāo)源和太陽(yáng)系內(nèi)的近地天體等快速運(yùn)動(dòng)目標(biāo).雖然這些天體位置與光度變化的時(shí)間尺度各不相同,且對(duì)于觀測(cè)提出的具體要求也各有不同,但是總體上都要求在各自的時(shí)間窗口內(nèi)保證觀測(cè)數(shù)據(jù)的連續(xù)獲取,因此大多使用多臺(tái)地理上分散分布的中小口徑望遠(yuǎn)鏡,通過(guò)自動(dòng)控制系統(tǒng)對(duì)天空進(jìn)行連續(xù)成像觀測(cè)[1?7].但這種觀測(cè)手段給時(shí)域天文的觀測(cè)和數(shù)據(jù)處理帶來(lái)了許多新的挑戰(zhàn)和問(wèn)題[8]:在時(shí)域天文觀測(cè)時(shí),除了近地快速移動(dòng)目標(biāo)為線源外,大部分感興趣的目標(biāo)屬于點(diǎn)源,在儀器端獲取的大量的天文數(shù)據(jù)中,會(huì)有部分?jǐn)?shù)據(jù)受到天空中云層的污染.在對(duì)這些數(shù)據(jù)中的天文目標(biāo)進(jìn)行探測(cè)和光度測(cè)量時(shí),云層往往會(huì)影響觀測(cè)精度,甚至?xí)耆蓴_觀測(cè)[9].如圖1所示為儀器端獲取的圖像數(shù)據(jù),其中(a)–(c)為沒(méi)有云的正常圖像,(d)–(f)為受云干擾的圖像.

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

由于云的影響,光透過(guò)云層強(qiáng)度會(huì)發(fā)生衰減,給測(cè)光帶來(lái)誤差;同時(shí)來(lái)自地面的光會(huì)從云區(qū)直接反射到成像系統(tǒng),使有云區(qū)比無(wú)云區(qū)更亮一些,嚴(yán)重影響對(duì)暗弱暫現(xiàn)源的提取.基于此,為了給測(cè)光和暗弱暫現(xiàn)源提取時(shí)提供參考,有必要首先對(duì)圖像中的云進(jìn)行提取并根據(jù)云的灰度變化建立指標(biāo)圖.指標(biāo)圖是一幅與原始數(shù)據(jù)大小相等的圖片,其每個(gè)像元的灰度信息可以反映出獲取的數(shù)據(jù)中的云的輪廓和云的灰度信息.但是如果圖像中沒(méi)有云,獲得的指標(biāo)圖對(duì)于觀測(cè)將變得毫無(wú)意義.因此,需要在提取云之前篩選出含云的圖像.由于時(shí)域天文觀測(cè)的數(shù)據(jù)量大,只依靠人工篩選耗時(shí)且繁重,所以有必要構(gòu)建云圖快速分類(lèi)系統(tǒng).之后根據(jù)分類(lèi)系統(tǒng)進(jìn)一步提取圖像中的云層,并建立指標(biāo)圖以便于科學(xué)研究[10?11].

但是,圖像的數(shù)據(jù)特性評(píng)價(jià)一直是一個(gè)非常困難的問(wèn)題,對(duì)于天文圖像尤甚.由于外界干擾因素很多,且不同的干擾因素對(duì)于不同天文觀測(cè)任務(wù)的影響也各不相同.因此,從成像的物理過(guò)程分析,通過(guò)確定幾個(gè)經(jīng)驗(yàn)參數(shù)及其分布范圍來(lái)直接判斷圖像中是否存在云非常困難.除了從物理過(guò)程直接分析外,另一種思路就是從數(shù)據(jù)特征出發(fā),將存在云的天文圖像進(jìn)行標(biāo)記,結(jié)合云圖像的形態(tài)特征,通過(guò)機(jī)器學(xué)習(xí)的方法建立符合數(shù)據(jù)特征的分類(lèi)器.

在星系分類(lèi)和光譜識(shí)別中[12?13],支持向量機(jī)(Support Vector Machine,SVM)是一種廣泛應(yīng)用的分類(lèi)器.該分類(lèi)器的基本策略是保證不同類(lèi)別的數(shù)據(jù)具有最大的分類(lèi)間隔.由于這類(lèi)求間隔最大化的問(wèn)題往往可以轉(zhuǎn)化為凸二次規(guī)劃問(wèn)題,因此與神經(jīng)網(wǎng)絡(luò)、隨機(jī)森林和決策樹(shù)等工具相比,SVM可以在數(shù)據(jù)量較少的情況下快速得到需要的分類(lèi)器[14?15].這一特性降低了數(shù)據(jù)積累的要求,同時(shí)減少了人工設(shè)置標(biāo)簽的工作量.此外,當(dāng)數(shù)據(jù)在低維空間不可分時(shí),通過(guò)核函數(shù)映射可以將分類(lèi)問(wèn)題轉(zhuǎn)化為求取高維空間分離超平面的問(wèn)題.本文將在第2章嘗試從云圖像的特征出發(fā),通過(guò)分析云圖像的特征,確定特征描述量構(gòu)建分類(lèi)系統(tǒng),并進(jìn)行了實(shí)驗(yàn)測(cè)試;第3章對(duì)分類(lèi)器獲得的有云的圖像進(jìn)行了指標(biāo)圖提取;最后對(duì)文章進(jìn)行了總結(jié).

2 天文云圖快速分類(lèi)系統(tǒng)

由于云是一種具有2維結(jié)構(gòu)的圖像,我們需要從它的特征出發(fā),利用一些參數(shù)對(duì)其進(jìn)行描述.一般來(lái)說(shuō),云層具有一定的延展特征和紋理特性,且云層的厚薄會(huì)導(dǎo)致灰度差異,為此,我們把圖像的灰度不一致度和紋理特征設(shè)定為天文圖像云層的評(píng)價(jià)指標(biāo),采用灰度不一致度和紋理特征共5個(gè)指標(biāo)共同對(duì)云進(jìn)行識(shí)別.

2.1 云圖特征提取

2.1.1 灰度不一致度

圖像的灰度不一致度是反映圖像背景灰度的一個(gè)重要指標(biāo),作為天文圖像的一個(gè)評(píng)價(jià)指標(biāo),圖像的灰度不一致度G的表達(dá)式如下:

式中,Ps和Pn分別為用局部圖像均方差的最大值與最小值.根據(jù)實(shí)際觀測(cè)時(shí)不同設(shè)備的工作情況可以得到典型的星像尺度,再根據(jù)其大小調(diào)節(jié)模板尺寸,本文選用9×9的模板遍歷全圖.

根據(jù)上述圖像的灰度不一致度的定義,我們可以發(fā)現(xiàn)該特征量所反映的是圖像在特征尺度內(nèi),灰度值整體波動(dòng)(背景)和局部特殊值(恒星或其他天文目標(biāo))的不一致度.圖像的灰度不一致度的高低對(duì)圖像質(zhì)量有著重要的參考價(jià)值,當(dāng)有云存在時(shí),圖像的灰度不一致度會(huì)發(fā)生劇烈變化.

2.1.2 紋理特征

在圖像灰度不一致度的基礎(chǔ)上,采用統(tǒng)計(jì)信息——紋理特征可以更好地對(duì)云進(jìn)一步描述.灰度共生矩陣能夠反映圖像灰度關(guān)于方向、相鄰間隔、變化幅度的綜合信息,所以其適用于對(duì)云圖像進(jìn)行紋理識(shí)別,進(jìn)而可以用于檢測(cè)天文圖像中是否有云的存在.滿(mǎn)足一定空間關(guān)系的灰度共生矩陣元素為:

式中,i和j分別為矩陣g的行數(shù)和列數(shù),M和N為矩陣g的行列數(shù),g為圖像f的灰度共生矩陣,#(x)表示集合中x元素的個(gè)數(shù).若(x1,z1)與(x2,z2)為矩陣f中的兩個(gè)點(diǎn),兩點(diǎn)之間距離為d,兩者與坐標(biāo)橫軸的夾角為q,則可以得到各種間距及角度的灰度共生矩陣g(i,j,d,q).我們這里研究的對(duì)象是天文圖像,主要的觀測(cè)目標(biāo)圖像一般都具有旋轉(zhuǎn)對(duì)稱(chēng)的特點(diǎn),這一特點(diǎn)導(dǎo)致在不同方向的灰度共生矩陣對(duì)圖像的紋理特征影響差異很小,所以這里選用d=1,q=0.考慮到望遠(yuǎn)鏡自動(dòng)觀測(cè)時(shí),獲取的天文圖像中云的特征,我們選用如下4個(gè)紋理特征對(duì)圖像進(jìn)行評(píng)價(jià):

(1)能量能量即每個(gè)元素的平方和,是對(duì)圖像平整度的衡量.如果g中的所有值均勻,則ASM的值較小;相反,在有云的情況下,部分區(qū)域有較大的值而其他區(qū)域值較小,則ASM的值較大,所以ASM可以用作對(duì)云圖的一個(gè)識(shí)別指標(biāo).

其中,Mg和Ng分別為矩陣g的行列數(shù).

(2)對(duì)比度對(duì)比度能夠反映圖像局部灰度變化的情況,云的存在會(huì)導(dǎo)致部分區(qū)域灰度變化很大,故有云圖像的CON值較大,所以CON可以作為云圖的另一個(gè)識(shí)別指標(biāo).

(3)逆差距逆差距反映了圖像紋理的同質(zhì)性,可以用來(lái)度量圖像紋理局部變化的多少.由于云往往具有一定的延展結(jié)構(gòu),其紋理特性變化會(huì)比較大,所以IDM可以作為檢測(cè)云圖的又一重要指標(biāo).

(4)熵熵是圖像所具有的信息量的度量,是對(duì)圖像隨機(jī)性的度量.當(dāng)相似觀測(cè)條件下的天文圖像存在云時(shí),熵值較不存在云時(shí)高.由于時(shí)域天文觀測(cè)圖像主要是點(diǎn)狀目標(biāo),其熵值小于同樣條件下包含云的圖像.因此ENT可以作為檢測(cè)云的最后一個(gè)重要指標(biāo).

2.2 云圖分類(lèi)器

2.2.1 SVM

SVM是由Vapnik等人提出的一種常用的分類(lèi)算法[16].最早被用于解決如下的二分類(lèi)問(wèn)題:

其中,w=(w1,w2,···,wn)為決定超平面方向的法向量,n為x的維度,x表示數(shù)據(jù)點(diǎn),b為超平面與原點(diǎn)之間的距離,SVM通過(guò)在高維數(shù)據(jù)集中尋找一個(gè)超平面來(lái)實(shí)現(xiàn)不同類(lèi)別之間的數(shù)據(jù)幾何間隔距離最大化.幾何間隔s定義為:

其中,y表示這些數(shù)據(jù)點(diǎn)代表的類(lèi)別,分別用1和–1表示,l為函數(shù)間隔,可以由下式求出:

則求幾何間隔距離最大化問(wèn)題可以轉(zhuǎn)化為下面的凸優(yōu)化問(wèn)題:

上述問(wèn)題可以使用二次規(guī)劃優(yōu)化包直接求解[17].但是,從上一小節(jié)的圖像特征描述模型討論中,我們發(fā)現(xiàn)各個(gè)評(píng)價(jià)指標(biāo)之間并不是獨(dú)立的,由于圖像的復(fù)雜性,這些指標(biāo)往往兩兩之間存在著關(guān)聯(lián).當(dāng)云層較薄時(shí),隨著云層面積增加,灰度整體波動(dòng)較大,ASM會(huì)加大,同時(shí)由于部分區(qū)域灰度變大,CON也會(huì)變化;但是,如果云層同時(shí)厚度增加,部分區(qū)域灰度不均勻性將降低,此時(shí)CON對(duì)應(yīng)值將會(huì)不變,甚至變小.由于這些量關(guān)聯(lián)性比較復(fù)雜,在可以滿(mǎn)足需求的情況下,如果再增加其他參量,關(guān)聯(lián)性可能會(huì)破壞特征空間的結(jié)構(gòu),導(dǎo)致SVM無(wú)法找到超平面,造成圖像分類(lèi)能力不佳.對(duì)于這一問(wèn)題,SVM可以將線性特征映射到更高維特征空間來(lái)實(shí)現(xiàn)圖像特征的可分[15].一個(gè)比較簡(jiǎn)單的辦法就是通過(guò)核函數(shù)的技巧,將上述分類(lèi)函數(shù)轉(zhuǎn)化為:

其中,K(xi,xj)為核函數(shù),ai為拉格朗日乘子,經(jīng)過(guò)轉(zhuǎn)化,凸優(yōu)化問(wèn)題成為了如下的形式:

2.2.2 SVM云圖分類(lèi)器構(gòu)建過(guò)程

我們首先使用線性核作為SVM核函數(shù),根據(jù)其對(duì)測(cè)試集的分類(lèi)性能進(jìn)行評(píng)價(jià).當(dāng)系統(tǒng)分類(lèi)性能較差時(shí)(準(zhǔn)確率小于90%),再進(jìn)一步考慮將SVM核轉(zhuǎn)換為非線性核.分類(lèi)系統(tǒng)流程如圖2所示.

從圖2可以看出,整個(gè)云圖像實(shí)時(shí)篩選系統(tǒng)包括如下步驟:

(1)讀取每一幅原始圖像的數(shù)據(jù),并人工根據(jù)原始數(shù)據(jù)是否有云添加標(biāo)簽;

(2)將讀取的數(shù)據(jù)分為兩部分,一部分用作訓(xùn)練集,其余部分用作測(cè)試集.根據(jù)經(jīng)驗(yàn),訓(xùn)練集數(shù)目應(yīng)該遠(yuǎn)大于特征維度,一般用100幅左右圖像可以達(dá)到訓(xùn)練的目的;

(3)分別對(duì)訓(xùn)練集和測(cè)試集的5個(gè)維度進(jìn)行特征提取,并對(duì)提取的特征歸一化處理.特征空間維度如表1所示,其中Original feature是原始的特征值,Normalized feature為規(guī)范化的特征值,我們選X1作為參考,將所有特征值映射到和X1同一個(gè)量級(jí),其中X2歸一化為10000X1,X3歸一化為L(zhǎng)n(X3),X4歸一化為100X4,X5保持原始值不變;

(4)用訓(xùn)練集提取的特征量訓(xùn)練SVM分類(lèi)器;

(5)用訓(xùn)練好的SVM分類(lèi)器對(duì)測(cè)試集進(jìn)行分類(lèi).X1–X5的歸一化是為了消除不同特征量級(jí)差距的影響,使SVM訓(xùn)練時(shí),保證數(shù)據(jù)在同一量級(jí).需要注意的是:歸一化是因數(shù)據(jù)而異的,也就是對(duì)于其他觀測(cè)需要根據(jù)數(shù)據(jù)特性設(shè)計(jì)其他歸一化方式,保證它們?cè)诹考?jí)上不會(huì)有差異[18].對(duì)圖1中的示例圖像,從左到右分別計(jì)算圖像的特征值,計(jì)算的特征值如表2所示.

圖2 系統(tǒng)流程圖Fig.2 System flowchart

表1 特征空間的維度Table 1 Dimensions of feature space

2.3 GPU加速

實(shí)際工作時(shí),該圖像分類(lèi)器將加載于控制計(jì)算機(jī)上對(duì)CCD采集的數(shù)據(jù)進(jìn)行分類(lèi).但是在G提取中采用了全圖遍歷,使得運(yùn)行的速度比較慢.處理大小為4096×4096的圖像,對(duì)G這一特征提取就需要耗時(shí)300秒·幅?1,因此我們需要對(duì)算法進(jìn)一步加速.

GPU又稱(chēng)圖像處理器,具有強(qiáng)大的浮點(diǎn)運(yùn)算和并行計(jì)算能力,與CPU相比更適合處理大量的并行數(shù)據(jù).本文提出的辦法在計(jì)算灰度不一致度時(shí)由于任意兩次運(yùn)算之間兩兩獨(dú)立,所以適合采用GPU進(jìn)行加速.望遠(yuǎn)鏡拍攝的圖片大小為M×N,遍歷采用的模板大小為m×n,則對(duì)單幅圖片需要運(yùn)算的次數(shù)為(M?m+1)×(N?n+1).

經(jīng)過(guò)測(cè)試,我們發(fā)現(xiàn)使用GPU加速后,對(duì)G這一特征的提取速度達(dá)到了0.43秒·幅?1,加速比接近700倍.

表2 正常圖像與含云圖像的特征值對(duì)比Table 2 The contrast of eigenvalues between normal images and images with cloud

2.4 分類(lèi)實(shí)驗(yàn)結(jié)果與分析

2.4.1 實(shí)驗(yàn)數(shù)據(jù)

本文實(shí)驗(yàn)數(shù)據(jù)選用的是南京大學(xué)時(shí)域天文臺(tái)(Time Domain Observatory,TIDO)18 cm望遠(yuǎn)鏡在紫金山天文臺(tái)盱眙觀測(cè)站試觀測(cè)時(shí)的數(shù)據(jù).望遠(yuǎn)鏡視場(chǎng)為12?,使用安道爾IkonXL系列CCD,像元數(shù)為4096×4096.我們從中選取665幅圖像作實(shí)驗(yàn)樣本,對(duì)這些數(shù)據(jù)進(jìn)行人工標(biāo)記,其中有云的圖片150幅,沒(méi)有云的515幅.對(duì)圖像有云、沒(méi)云分別標(biāo)注為類(lèi)1、類(lèi)2.

2.4.2 實(shí)驗(yàn)過(guò)程與結(jié)果

我們分別從類(lèi)1、類(lèi)2隨機(jī)選取50%作為訓(xùn)練集,剩余的作為測(cè)試集數(shù)據(jù).然后按照?qǐng)D2中所示的流程開(kāi)始進(jìn)行實(shí)驗(yàn)測(cè)試,最終獲得的結(jié)果用圖3的混淆矩陣表示出來(lái).

圖3 混淆矩陣Fig.3 The confusion matrix

從圖3中可以看出:在332幅測(cè)試圖像中,類(lèi)1共有75幅,在實(shí)際測(cè)試中有1張圖像被錯(cuò)分為類(lèi)2,其錯(cuò)誤率為1.3%;類(lèi)2共有257幅,在實(shí)際測(cè)試中有3張被誤分為類(lèi)1,其錯(cuò)誤率為1.2%.類(lèi)1、類(lèi)2的總體錯(cuò)誤率為1.2%.這一分類(lèi)準(zhǔn)確率基本達(dá)到了我們的要求,為下一步提取云提供了參考.

3 云的提取

為了保證觀測(cè)的空間和時(shí)間連續(xù)性,雖然含云的圖像質(zhì)量下降,但通過(guò)參考指標(biāo)圖,去除掉觀測(cè)數(shù)據(jù)中含云的部分,剩余的數(shù)據(jù)仍然可以利用.圖4為我們?cè)O(shè)計(jì)的分類(lèi)算法識(shí)別出的有云的圖像,下面以圖4為例介紹云的指標(biāo)圖的提取過(guò)程.

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

由于云是整體具有一定輪廓和紋理的圖像,從圖像信息的頻率域分析,其信息主要分布于中低頻;從圖像像元的灰度分布分析,當(dāng)云的灰度值比較高時(shí)(比較厚),包含云的像元會(huì)在直方圖上產(chǎn)生一個(gè)峰,當(dāng)云的灰度值比較低時(shí)(比較薄),包含云的像元會(huì)擴(kuò)展背景像元形成的峰的寬度.為此,我們?cè)O(shè)計(jì)了如下過(guò)程對(duì)圖像中含云的部分進(jìn)行提取:

(1)提取圖像的灰度直方圖,對(duì)圖像進(jìn)行高斯多峰擬合,并求擬合函數(shù)的拐點(diǎn).由于圖像的灰度值主要分布在圖像的背景區(qū)域,圖像灰度直方圖近似高斯分布,故采用如下多高斯函數(shù)進(jìn)行擬合:

其中,a1、b1、c1、a2、b2、c2、a3、b3、c3為擬合函數(shù)的系數(shù), 求擬合函數(shù)系數(shù),對(duì)h求2階導(dǎo)數(shù),與x軸的交點(diǎn)即為函數(shù)的拐點(diǎn),分別取兩個(gè)拐點(diǎn)值為背景灰度值的上下閾值,圖5為灰度直方圖的多高斯擬合曲線h.

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

(2)分別取兩個(gè)拐點(diǎn)值為背景灰度值的上下閾值,將這兩個(gè)閾值之間的部分6等分,依次用不同的閾值繪制圖像等高線并刪除面積小于500的區(qū)域,用該閾值替代這一梯度的灰度值,所得的灰度梯度圖如圖6所示.

(3)對(duì)原始圖像進(jìn)行濾波,由于圖像內(nèi)背景和云是低頻信息,而圖像中的星像為高頻信息,通過(guò)使用低通濾波器,可以降低個(gè)別亮星的干擾.為此我們選擇高斯濾波,通過(guò)SExtractor提取星的位置信息,預(yù)估圖像中星像的最大直徑,結(jié)合經(jīng)驗(yàn)選擇合適大小的核對(duì)圖像進(jìn)行濾波,我們這里核大小選擇為50.

(4)對(duì)濾波后的圖像,按照(2)中所求的梯度,分別在該梯度進(jìn)行如下操作:

所求的云的指標(biāo)圖如圖7所示.從圖中可以看到,大部分云的輪廓可以提取出來(lái),能夠?yàn)闇y(cè)光和后續(xù)暗弱目標(biāo)提取的參考.但是亮星的影響比較難去除,是影響這個(gè)辦法性能的主要因素.圖8是更多的含云圖像處理結(jié)果,其中(a)–(c)為含云圖像,(d)–(f)為云的指標(biāo)圖.

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

圖7 云的指標(biāo)圖Fig.7 The index figure of clouds

圖8 含云圖像及其處理結(jié)果.(a)–(c)含云圖像;(d)–(f)云的指標(biāo)圖Fig.8 Images with colud and their process result.(a)–(c)images with colud;(d)–(f)the index figure of clouds

4 總結(jié)

本文根據(jù)時(shí)域天文觀測(cè)數(shù)據(jù)處理需求,針對(duì)天文圖像中云的特點(diǎn)提出了融合灰度不一致度與灰度紋理的圖像特征,并據(jù)此構(gòu)建了基于SVM的實(shí)時(shí)云圖篩選系統(tǒng),在篩選獲得的圖像基礎(chǔ)上對(duì)圖片中云的輪廓進(jìn)行了提取.該系統(tǒng)的分類(lèi)準(zhǔn)確度可以達(dá)到98%以上,準(zhǔn)確率高,魯棒性強(qiáng),同時(shí)也彌補(bǔ)了人工篩選中自動(dòng)化程度低的不足.從實(shí)際應(yīng)用角度,系統(tǒng)中采用了GPU技術(shù),滿(mǎn)足高速處理的需求.經(jīng)過(guò)處理后獲得的云的指標(biāo)圖可以大致反映云的輪廓,但是部分結(jié)果會(huì)受到亮星的影響,這一問(wèn)題有待于進(jìn)一步研究.

致謝感謝審稿人提出的建議和紫金山天文臺(tái)孫榮煜對(duì)本研究的大力支持.