王利文 賈 鵬? 蔡冬梅 劉慧根
(1太原理工大學(xué)物理與光電工程學(xué)院太原030024)
(2南京大學(xué)天文與空間科學(xué)學(xué)院南京210034)
時(shí)域天文學(xué)是目前非?;钴S的一個(gè)研究領(lǐng)域,其研究對象主要包括超新星、變星及系外行星等光度連續(xù)變化的目標(biāo)源和太陽系內(nèi)的近地天體等快速運(yùn)動(dòng)目標(biāo).雖然這些天體位置與光度變化的時(shí)間尺度各不相同,且對于觀測提出的具體要求也各有不同,但是總體上都要求在各自的時(shí)間窗口內(nèi)保證觀測數(shù)據(jù)的連續(xù)獲取,因此大多使用多臺地理上分散分布的中小口徑望遠(yuǎn)鏡,通過自動(dòng)控制系統(tǒng)對天空進(jìn)行連續(xù)成像觀測[1?7].但這種觀測手段給時(shí)域天文的觀測和數(shù)據(jù)處理帶來了許多新的挑戰(zhàn)和問題[8]:在時(shí)域天文觀測時(shí),除了近地快速移動(dòng)目標(biāo)為線源外,大部分感興趣的目標(biāo)屬于點(diǎn)源,在儀器端獲取的大量的天文數(shù)據(jù)中,會(huì)有部分?jǐn)?shù)據(jù)受到天空中云層的污染.在對這些數(shù)據(jù)中的天文目標(biāo)進(jìn)行探測和光度測量時(shí),云層往往會(huì)影響觀測精度,甚至?xí)耆蓴_觀測[9].如圖1所示為儀器端獲取的圖像數(shù)據(jù),其中(a)–(c)為沒有云的正常圖像,(d)–(f)為受云干擾的圖像.
圖1 圖像示例.(a)–(c)正常圖像;(d)–(f)含云圖像Fig.1 Image examples.(a)–(c)normal images;(d)–(f)images with cloud
由于云的影響,光透過云層強(qiáng)度會(huì)發(fā)生衰減,給測光帶來誤差;同時(shí)來自地面的光會(huì)從云區(qū)直接反射到成像系統(tǒng),使有云區(qū)比無云區(qū)更亮一些,嚴(yán)重影響對暗弱暫現(xiàn)源的提取.基于此,為了給測光和暗弱暫現(xiàn)源提取時(shí)提供參考,有必要首先對圖像中的云進(jìn)行提取并根據(jù)云的灰度變化建立指標(biāo)圖.指標(biāo)圖是一幅與原始數(shù)據(jù)大小相等的圖片,其每個(gè)像元的灰度信息可以反映出獲取的數(shù)據(jù)中的云的輪廓和云的灰度信息.但是如果圖像中沒有云,獲得的指標(biāo)圖對于觀測將變得毫無意義.因此,需要在提取云之前篩選出含云的圖像.由于時(shí)域天文觀測的數(shù)據(jù)量大,只依靠人工篩選耗時(shí)且繁重,所以有必要構(gòu)建云圖快速分類系統(tǒng).之后根據(jù)分類系統(tǒng)進(jìn)一步提取圖像中的云層,并建立指標(biāo)圖以便于科學(xué)研究[10?11].
但是,圖像的數(shù)據(jù)特性評價(jià)一直是一個(gè)非常困難的問題,對于天文圖像尤甚.由于外界干擾因素很多,且不同的干擾因素對于不同天文觀測任務(wù)的影響也各不相同.因此,從成像的物理過程分析,通過確定幾個(gè)經(jīng)驗(yàn)參數(shù)及其分布范圍來直接判斷圖像中是否存在云非常困難.除了從物理過程直接分析外,另一種思路就是從數(shù)據(jù)特征出發(fā),將存在云的天文圖像進(jìn)行標(biāo)記,結(jié)合云圖像的形態(tài)特征,通過機(jī)器學(xué)習(xí)的方法建立符合數(shù)據(jù)特征的分類器.
在星系分類和光譜識別中[12?13],支持向量機(jī)(Support Vector Machine,SVM)是一種廣泛應(yīng)用的分類器.該分類器的基本策略是保證不同類別的數(shù)據(jù)具有最大的分類間隔.由于這類求間隔最大化的問題往往可以轉(zhuǎn)化為凸二次規(guī)劃問題,因此與神經(jīng)網(wǎng)絡(luò)、隨機(jī)森林和決策樹等工具相比,SVM可以在數(shù)據(jù)量較少的情況下快速得到需要的分類器[14?15].這一特性降低了數(shù)據(jù)積累的要求,同時(shí)減少了人工設(shè)置標(biāo)簽的工作量.此外,當(dāng)數(shù)據(jù)在低維空間不可分時(shí),通過核函數(shù)映射可以將分類問題轉(zhuǎn)化為求取高維空間分離超平面的問題.本文將在第2章嘗試從云圖像的特征出發(fā),通過分析云圖像的特征,確定特征描述量構(gòu)建分類系統(tǒng),并進(jìn)行了實(shí)驗(yàn)測試;第3章對分類器獲得的有云的圖像進(jìn)行了指標(biāo)圖提取;最后對文章進(jìn)行了總結(jié).
由于云是一種具有2維結(jié)構(gòu)的圖像,我們需要從它的特征出發(fā),利用一些參數(shù)對其進(jìn)行描述.一般來說,云層具有一定的延展特征和紋理特性,且云層的厚薄會(huì)導(dǎo)致灰度差異,為此,我們把圖像的灰度不一致度和紋理特征設(shè)定為天文圖像云層的評價(jià)指標(biāo),采用灰度不一致度和紋理特征共5個(gè)指標(biāo)共同對云進(jìn)行識別.
2.1.1 灰度不一致度
圖像的灰度不一致度是反映圖像背景灰度的一個(gè)重要指標(biāo),作為天文圖像的一個(gè)評價(jià)指標(biāo),圖像的灰度不一致度G的表達(dá)式如下:
式中,Ps和Pn分別為用局部圖像均方差的最大值與最小值.根據(jù)實(shí)際觀測時(shí)不同設(shè)備的工作情況可以得到典型的星像尺度,再根據(jù)其大小調(diào)節(jié)模板尺寸,本文選用9×9的模板遍歷全圖.
根據(jù)上述圖像的灰度不一致度的定義,我們可以發(fā)現(xiàn)該特征量所反映的是圖像在特征尺度內(nèi),灰度值整體波動(dòng)(背景)和局部特殊值(恒星或其他天文目標(biāo))的不一致度.圖像的灰度不一致度的高低對圖像質(zhì)量有著重要的參考價(jià)值,當(dāng)有云存在時(shí),圖像的灰度不一致度會(huì)發(fā)生劇烈變化.
2.1.2 紋理特征
在圖像灰度不一致度的基礎(chǔ)上,采用統(tǒng)計(jì)信息——紋理特征可以更好地對云進(jìn)一步描述.灰度共生矩陣能夠反映圖像灰度關(guān)于方向、相鄰間隔、變化幅度的綜合信息,所以其適用于對云圖像進(jìn)行紋理識別,進(jì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).我們這里研究的對象是天文圖像,主要的觀測目標(biāo)圖像一般都具有旋轉(zhuǎn)對稱的特點(diǎn),這一特點(diǎn)導(dǎo)致在不同方向的灰度共生矩陣對圖像的紋理特征影響差異很小,所以這里選用d=1,q=0.考慮到望遠(yuǎn)鏡自動(dòng)觀測時(shí),獲取的天文圖像中云的特征,我們選用如下4個(gè)紋理特征對圖像進(jìn)行評價(jià):
(1)能量能量即每個(gè)元素的平方和,是對圖像平整度的衡量.如果g中的所有值均勻,則ASM的值較小;相反,在有云的情況下,部分區(qū)域有較大的值而其他區(qū)域值較小,則ASM的值較大,所以ASM可以用作對云圖的一個(gè)識別指標(biāo).
其中,Mg和Ng分別為矩陣g的行列數(shù).
(2)對比度對比度能夠反映圖像局部灰度變化的情況,云的存在會(huì)導(dǎo)致部分區(qū)域灰度變化很大,故有云圖像的CON值較大,所以CON可以作為云圖的另一個(gè)識別指標(biāo).
(3)逆差距逆差距反映了圖像紋理的同質(zhì)性,可以用來度量圖像紋理局部變化的多少.由于云往往具有一定的延展結(jié)構(gòu),其紋理特性變化會(huì)比較大,所以IDM可以作為檢測云圖的又一重要指標(biāo).
(4)熵熵是圖像所具有的信息量的度量,是對圖像隨機(jī)性的度量.當(dāng)相似觀測條件下的天文圖像存在云時(shí),熵值較不存在云時(shí)高.由于時(shí)域天文觀測圖像主要是點(diǎn)狀目標(biāo),其熵值小于同樣條件下包含云的圖像.因此ENT可以作為檢測云的最后一個(gè)重要指標(biāo).
2.2.1 SVM
SVM是由Vapnik等人提出的一種常用的分類算法[16].最早被用于解決如下的二分類問題:
其中,w=(w1,w2,···,wn)為決定超平面方向的法向量,n為x的維度,x表示數(shù)據(jù)點(diǎn),b為超平面與原點(diǎn)之間的距離,SVM通過在高維數(shù)據(jù)集中尋找一個(gè)超平面來實(shí)現(xiàn)不同類別之間的數(shù)據(jù)幾何間隔距離最大化.幾何間隔s定義為:
其中,y表示這些數(shù)據(jù)點(diǎn)代表的類別,分別用1和–1表示,l為函數(shù)間隔,可以由下式求出:
則求幾何間隔距離最大化問題可以轉(zhuǎn)化為下面的凸優(yōu)化問題:
上述問題可以使用二次規(guī)劃優(yōu)化包直接求解[17].但是,從上一小節(jié)的圖像特征描述模型討論中,我們發(fā)現(xiàn)各個(gè)評價(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對應(yīng)值將會(huì)不變,甚至變小.由于這些量關(guān)聯(lián)性比較復(fù)雜,在可以滿足需求的情況下,如果再增加其他參量,關(guān)聯(lián)性可能會(huì)破壞特征空間的結(jié)構(gòu),導(dǎo)致SVM無法找到超平面,造成圖像分類能力不佳.對于這一問題,SVM可以將線性特征映射到更高維特征空間來實(shí)現(xiàn)圖像特征的可分[15].一個(gè)比較簡單的辦法就是通過核函數(shù)的技巧,將上述分類函數(shù)轉(zhuǎn)化為:
其中,K(xi,xj)為核函數(shù),ai為拉格朗日乘子,經(jīng)過轉(zhuǎn)化,凸優(yōu)化問題成為了如下的形式:
2.2.2 SVM云圖分類器構(gòu)建過程
我們首先使用線性核作為SVM核函數(shù),根據(jù)其對測試集的分類性能進(jìn)行評價(jià).當(dāng)系統(tǒng)分類性能較差時(shí)(準(zhǔn)確率小于90%),再進(jìn)一步考慮將SVM核轉(zhuǎn)換為非線性核.分類系統(tǒng)流程如圖2所示.
從圖2可以看出,整個(gè)云圖像實(shí)時(shí)篩選系統(tǒng)包括如下步驟:
(1)讀取每一幅原始圖像的數(shù)據(jù),并人工根據(jù)原始數(shù)據(jù)是否有云添加標(biāo)簽;
(2)將讀取的數(shù)據(jù)分為兩部分,一部分用作訓(xùn)練集,其余部分用作測試集.根據(jù)經(jīng)驗(yàn),訓(xùn)練集數(shù)目應(yīng)該遠(yuǎn)大于特征維度,一般用100幅左右圖像可以達(dá)到訓(xùn)練的目的;
(3)分別對訓(xùn)練集和測試集的5個(gè)維度進(jìn)行特征提取,并對提取的特征歸一化處理.特征空間維度如表1所示,其中Original feature是原始的特征值,Normalized feature為規(guī)范化的特征值,我們選X1作為參考,將所有特征值映射到和X1同一個(gè)量級,其中X2歸一化為10000X1,X3歸一化為Ln(X3),X4歸一化為100X4,X5保持原始值不變;
(4)用訓(xùn)練集提取的特征量訓(xùn)練SVM分類器;
(5)用訓(xùn)練好的SVM分類器對測試集進(jìn)行分類.X1–X5的歸一化是為了消除不同特征量級差距的影響,使SVM訓(xùn)練時(shí),保證數(shù)據(jù)在同一量級.需要注意的是:歸一化是因數(shù)據(jù)而異的,也就是對于其他觀測需要根據(jù)數(shù)據(jù)特性設(shè)計(jì)其他歸一化方式,保證它們在量級上不會(huì)有差異[18].對圖1中的示例圖像,從左到右分別計(jì)算圖像的特征值,計(jì)算的特征值如表2所示.
圖2 系統(tǒng)流程圖Fig.2 System flowchart
表1 特征空間的維度Table 1 Dimensions of feature space
實(shí)際工作時(shí),該圖像分類器將加載于控制計(jì)算機(jī)上對CCD采集的數(shù)據(jù)進(jìn)行分類.但是在G提取中采用了全圖遍歷,使得運(yùn)行的速度比較慢.處理大小為4096×4096的圖像,對G這一特征提取就需要耗時(shí)300秒·幅?1,因此我們需要對算法進(jìn)一步加速.
GPU又稱圖像處理器,具有強(qiáng)大的浮點(diǎn)運(yùn)算和并行計(jì)算能力,與CPU相比更適合處理大量的并行數(shù)據(jù).本文提出的辦法在計(jì)算灰度不一致度時(shí)由于任意兩次運(yùn)算之間兩兩獨(dú)立,所以適合采用GPU進(jìn)行加速.望遠(yuǎn)鏡拍攝的圖片大小為M×N,遍歷采用的模板大小為m×n,則對單幅圖片需要運(yùn)算的次數(shù)為(M?m+1)×(N?n+1).
經(jīng)過測試,我們發(fā)現(xiàn)使用GPU加速后,對G這一特征的提取速度達(dá)到了0.43秒·幅?1,加速比接近700倍.
表2 正常圖像與含云圖像的特征值對比Table 2 The contrast of eigenvalues between normal images and images with cloud
2.4.1 實(shí)驗(yàn)數(shù)據(jù)
本文實(shí)驗(yàn)數(shù)據(jù)選用的是南京大學(xué)時(shí)域天文臺(Time Domain Observatory,TIDO)18 cm望遠(yuǎn)鏡在紫金山天文臺盱眙觀測站試觀測時(shí)的數(shù)據(jù).望遠(yuǎn)鏡視場為12?,使用安道爾IkonXL系列CCD,像元數(shù)為4096×4096.我們從中選取665幅圖像作實(shí)驗(yàn)樣本,對這些數(shù)據(jù)進(jìn)行人工標(biāo)記,其中有云的圖片150幅,沒有云的515幅.對圖像有云、沒云分別標(biāo)注為類1、類2.
2.4.2 實(shí)驗(yàn)過程與結(jié)果
我們分別從類1、類2隨機(jī)選取50%作為訓(xùn)練集,剩余的作為測試集數(shù)據(jù).然后按照圖2中所示的流程開始進(jìn)行實(shí)驗(yàn)測試,最終獲得的結(jié)果用圖3的混淆矩陣表示出來.
圖3 混淆矩陣Fig.3 The confusion matrix
從圖3中可以看出:在332幅測試圖像中,類1共有75幅,在實(shí)際測試中有1張圖像被錯(cuò)分為類2,其錯(cuò)誤率為1.3%;類2共有257幅,在實(shí)際測試中有3張被誤分為類1,其錯(cuò)誤率為1.2%.類1、類2的總體錯(cuò)誤率為1.2%.這一分類準(zhǔn)確率基本達(dá)到了我們的要求,為下一步提取云提供了參考.
為了保證觀測的空間和時(shí)間連續(xù)性,雖然含云的圖像質(zhì)量下降,但通過參考指標(biāo)圖,去除掉觀測數(shù)據(jù)中含云的部分,剩余的數(shù)據(jù)仍然可以利用.圖4為我們設(shè)計(jì)的分類算法識別出的有云的圖像,下面以圖4為例介紹云的指標(biāo)圖的提取過程.
圖4 有云圖像Fig.4 Image with clouds
由于云是整體具有一定輪廓和紋理的圖像,從圖像信息的頻率域分析,其信息主要分布于中低頻;從圖像像元的灰度分布分析,當(dāng)云的灰度值比較高時(shí)(比較厚),包含云的像元會(huì)在直方圖上產(chǎn)生一個(gè)峰,當(dāng)云的灰度值比較低時(shí)(比較薄),包含云的像元會(huì)擴(kuò)展背景像元形成的峰的寬度.為此,我們設(shè)計(jì)了如下過程對圖像中含云的部分進(jìn)行提取:
(1)提取圖像的灰度直方圖,對圖像進(jìn)行高斯多峰擬合,并求擬合函數(shù)的拐點(diǎn).由于圖像的灰度值主要分布在圖像的背景區(qū)域,圖像灰度直方圖近似高斯分布,故采用如下多高斯函數(shù)進(jìn)行擬合:
其中,a1、b1、c1、a2、b2、c2、a3、b3、c3為擬合函數(shù)的系數(shù), 求擬合函數(shù)系數(shù),對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)對原始圖像進(jìn)行濾波,由于圖像內(nèi)背景和云是低頻信息,而圖像中的星像為高頻信息,通過使用低通濾波器,可以降低個(gè)別亮星的干擾.為此我們選擇高斯濾波,通過SExtractor提取星的位置信息,預(yù)估圖像中星像的最大直徑,結(jié)合經(jīng)驗(yàn)選擇合適大小的核對圖像進(jìn)行濾波,我們這里核大小選擇為50.
(4)對濾波后的圖像,按照(2)中所求的梯度,分別在該梯度進(jìn)行如下操作:
所求的云的指標(biāo)圖如圖7所示.從圖中可以看到,大部分云的輪廓可以提取出來,能夠?yàn)闇y光和后續(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
本文根據(jù)時(shí)域天文觀測數(shù)據(jù)處理需求,針對天文圖像中云的特點(diǎn)提出了融合灰度不一致度與灰度紋理的圖像特征,并據(jù)此構(gòu)建了基于SVM的實(shí)時(shí)云圖篩選系統(tǒng),在篩選獲得的圖像基礎(chǔ)上對圖片中云的輪廓進(jìn)行了提取.該系統(tǒng)的分類準(zhǔn)確度可以達(dá)到98%以上,準(zhǔn)確率高,魯棒性強(qiáng),同時(shí)也彌補(bǔ)了人工篩選中自動(dòng)化程度低的不足.從實(shí)際應(yīng)用角度,系統(tǒng)中采用了GPU技術(shù),滿足高速處理的需求.經(jīng)過處理后獲得的云的指標(biāo)圖可以大致反映云的輪廓,但是部分結(jié)果會(huì)受到亮星的影響,這一問題有待于進(jìn)一步研究.
致謝感謝審稿人提出的建議和紫金山天文臺孫榮煜對本研究的大力支持.