基于NPSVM模型的超新星識別方法?

王精東 陳星星 梁吳穎 黃澤鋒 全少武 沈錦泫

石雄輝1?

(1 廣東海洋大學(xué)數(shù)學(xué)與計(jì)算機(jī)學(xué)院湛江524088)(2 北京師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院北京100875)

(3 湖南科技大學(xué)材料科學(xué)與工程學(xué)院湘潭411201)

1 引言

遙遠(yuǎn)星系內(nèi)的恒星在一般情況下亮度較低, 無法觀測, 而其爆發(fā)成為超新星時(shí)則突然進(jìn)入了可觀測的范圍, 此現(xiàn)象作為恒星演化的一個(gè)重要過程, 對研究恒星與高能天體物理具有重要影響[1]. 由于超新星的出現(xiàn)在時(shí)間與空間上均具有隨機(jī)性, 使得超新星觀測大多采用大規(guī)模巡天以獲得大量對象的樣本數(shù)據(jù)進(jìn)行分析篩選.

超新星搜尋主要有可見光圖像搜尋、測光、光譜認(rèn)證3種方法[1–2], 其中, 光譜數(shù)據(jù)相對可見光圖像, 具有更好的準(zhǔn)確性, 且不用像可見光圖像、測光一樣需獲取多日的數(shù)據(jù)進(jìn)行比對校驗(yàn), 是超新星搜尋的重要手段.

目前在運(yùn)行的大型超新星巡天項(xiàng)目中, 國外有ZTF (the Zwicky Transient Facility)、Pan-STARRS (Panoramic Survey Telescope And Rapid Response System)、Gaia (Global Astrometric Interferometer for Astrophysics)等[3–5], 國內(nèi)有TNTS (the THU-NAOC Transient Survey)、PTSS(the PMO-Tsinghua Supernova Survey)、AST3(Antarctic Survey Telescope)等項(xiàng)目[6–8]. 這些大型巡天計(jì)劃往往產(chǎn)生超大規(guī)模的數(shù)據(jù),現(xiàn)代天文觀測對自動化分析的需求日益明顯.

光譜數(shù)據(jù)往往由于維數(shù)較多, 從而需要降維以方便計(jì)算, 常用的方法有主成分分析(Principal Component Analysis, PCA), 比如, 劉真祥等人直接采用PCA降維[9]; 屠良平等人基于PCA實(shí)現(xiàn)光譜分解并計(jì)算重構(gòu)系數(shù)作為新特征[10]. 在識別方面, 可以采用深度信念網(wǎng)絡(luò), 劉真祥等人采用該方法在激變變星光譜分類問題上獲得了97.4%的正確率[9]; 屠良平等人則采用基于局部孤立性因子(Local Outlier Factor, LOF)的離群點(diǎn)檢測算法將超新星候選數(shù)降低到了實(shí)驗(yàn)總數(shù)的1%, 該實(shí)驗(yàn)每次隨機(jī)將1條超新星光譜混入到5054條星系光譜中, 發(fā)現(xiàn)每次實(shí)驗(yàn)中超新星的LOF都排在實(shí)驗(yàn)樣本集合的前1%[11].此外還可以使用支持向量機(jī)(Support Vector Machine, SVM)或人工神經(jīng)網(wǎng)絡(luò)(Artificial Neural Network, ANN)等方法[12–13].

經(jīng)研究發(fā)現(xiàn), 基于光譜分解的特征提取方法存在冗余特征的問題, 并提出了解決方案: 首先, 用光譜分解方法提取特征進(jìn)行初篩, 然后計(jì)算所得特征系數(shù)在兩個(gè)類別上的相似度, 并做進(jìn)一步篩選; SVM在星系-超新星二分類問題中得到了99.502%的識別正確率. 考慮進(jìn)一步降低漏判率, 結(jié)合Neyman-Pearson決策, 本文提出了NPSVM(Neyman-Pearson SVM)分類模型, 以適應(yīng)不同條件與場景.

2 數(shù)據(jù)收集

超新星往往是從星系對象中篩選出來. 所以只需要收集超新星與星系兩類的數(shù)據(jù)進(jìn)行分析即可.

我們從WISeREP (the Weizmann Interactive Supernova data REPository)1https://wiserep.weizmann.ac.il上收集了SUSPECT (SUpernova SPECTrum), SNfactory (the Nearby Supernova Factory),SDSS-SNe (the Sloan Digital Sky Survey-II Supernova Survey)觀測數(shù)據(jù)光譜[14], 使用了CfA (the Center for Astrophysics)2https://www.cfa.harvard.edu/supernova/cfaspecsnIa 20120322.tar.gz超新星計(jì)劃[15]的一部分與Nugent等人構(gòu)造的超新星模板3https://c3.lbl.gov/nugent/nugent templates.html[16]. 本文與文獻(xiàn)[10–11]中僅采用Ia型超新星模板不同, 希望能夠?qū)⒛Ｐ屯茝V到其他類型的超新星光譜識別上, 所以我們也應(yīng)用了Nugent等人的其他類型超新星模板,由于這些模板是由超新星同一時(shí)間的多個(gè)觀測來源的光譜取平均得到, 并在時(shí)間上進(jìn)行了插值[16], 更具有代表性, 因此將用于構(gòu)造超新星的PCA模板, 而剩下的超新星實(shí)測光譜將用于分類模型的訓(xùn)練和測試.

星系數(shù)據(jù)上, 我們使用了SDSS-dr16得到的星系觀測數(shù)據(jù)[17]. 這些星系光譜將既用于PCA模板庫的構(gòu)造又作為實(shí)驗(yàn)數(shù)據(jù)集的星系部分.

為了將數(shù)據(jù)統(tǒng)一在一個(gè)區(qū)間上進(jìn)行分析, 本文篩選出靜止波長范圍在3809–7000 ?A的光譜, 并獲得了3427條超新星實(shí)測光譜和2193條星系光譜.

表1、2分別給出了超新星和星系實(shí)測光譜的類別數(shù)量和比例, 未給出子類別的記為untyped, 其中AGN (Active Galactic Nucleus)為活動星系核. 圖1給出了超新星實(shí)測光譜相對光極大和數(shù)據(jù)集信噪比的主要分布, 圖例中SNe表示超新星, Gal表示星系.

表1 超新星(SN)類型統(tǒng)計(jì)Table 1 Statistics of supernova types

表2 星系類型統(tǒng)計(jì)Table 2 Statistics of galaxy types

3 數(shù)據(jù)預(yù)處理

超新星形成機(jī)制包括白矮星吸積伴星產(chǎn)生的熱核爆炸、大質(zhì)量恒星的核塌縮爆炸,兩者均多發(fā)生在星系內(nèi). 所以超新星的光譜S往往會與宿主星系的光譜G發(fā)生疊加, 同時(shí), 伴隨有測量、大氣條件、天光背景等引起的隨機(jī)噪聲ε, 于是觀測得到的光譜M可以寫成如下形式:

光譜作為對遙遠(yuǎn)天體的觀測手段之一, 除了噪聲影響外, 主要還存在以下影響因素:

(1)遙遠(yuǎn)天體的光譜紅移;

(2)天體距離、超新星本身屬性、所處演化階段等對譜線強(qiáng)度的影響;

(3)不同望遠(yuǎn)鏡數(shù)據(jù)處理方法的差異.

圖1 超新星相對光極大(a)和數(shù)據(jù)集信噪比(b)的分布直方圖. (a)中收集自CfA的光譜(紅色曲線)給出的是相對B波段光極大的天數(shù), 收集自WISeREP的光譜(藍(lán)色曲線)給出的是相對于V波段光極大的天數(shù). (b)中紅色曲線為超新星光譜的信噪比分布, 藍(lán)色曲線為星系光譜的信噪比分布.Fig.1 Histogram of supernova phases (a) and signal-to-noise ratio of dataset (b). In panel (a), the spectra collected from CfA give the number of days relative to B-band maximum (red curve), while the spectra collected from WISeREP give the number of days relative to V-band maximum (blue curve). In panel (b), the red curve is the signal-to-noise ratio distribution of supernova spectra, and the blue curve is the signal-to-noise ratio distribution of galactic spectra.

因此本文對觀測數(shù)據(jù)進(jìn)行了以下預(yù)處理以減小上述因素對算法效果的影響:

(1)退紅移: 對于紅移為z、觀測波長為λobs的天體光譜, 計(jì)算其靜止波長:

(2)截?cái)嗯c插值: 將所有樣本截?cái)嘀领o止波長覆蓋范圍的交集3809–7000 ?A并進(jìn)行線性插值;

(3)去除強(qiáng)窄線: 由于超新星的光譜是沒有窄線的[10], 為了減少窄線的干擾, 本文基于3σ原則(將到均值的距離大于3倍標(biāo)準(zhǔn)差的點(diǎn)視為異常點(diǎn))來檢測窄線, 對窄線進(jìn)行過濾并替換成均值, 循環(huán)多次以保證效果;

(4)去噪: 使用小波變換方法進(jìn)行去噪;

(5)標(biāo)準(zhǔn)化: 對于均值為μi、標(biāo)準(zhǔn)差為σi的第i條觀測數(shù)據(jù)Mobs,i, 標(biāo)準(zhǔn)化后記為向量Mi=(Mi1,Mi2,··· ,Miq). 其第j個(gè)分量:

其中,Mobs,ij為觀測數(shù)據(jù)向量Mobs,i的第j個(gè)分量,N為觀測樣本的數(shù)量,q為光譜數(shù)據(jù)向量的維數(shù). 圖2展示了數(shù)據(jù)預(yù)處理的效果.

圖2 在數(shù)據(jù)預(yù)處理中, 一條光譜的原始數(shù)據(jù)將被退紅移、截?cái)?、插值、去譜線、去噪、標(biāo)準(zhǔn)化.Fig.2 In the data preprocessing, the raw data of a spectrum will be deredshifted, truncated,interpolated, spectral lines removed, denoised and normalised.

4 特征提取

以下設(shè)行向量

為第i條經(jīng)過預(yù)處理的混合光譜, 其中,Gi與Si分別表示組成該混合光譜的星系光譜和超新星光譜,εi表示觀測過程中的隨機(jī)噪聲.

4.1 光譜分解

PCA旨在通過線性變換得到新的變量, 即主成分, 并對主成分按照方差貢獻(xiàn)率進(jìn)行排序. 這些主成分滿足: 各個(gè)主成分Yj的方差達(dá)到最大且主成分之間協(xié)方差為0,即第j個(gè)主成分Yj的方差Var(Yj) = max,Yj與Yj′的協(xié)方差Cov(Yj,Yj′) = 0, j >1, j′=1,2,...,j ?1. 對于第i條光譜Mi, 其第j個(gè)主成分其中,lj稱為主成分方向.

對超新星與星系模板庫分別運(yùn)用PCA, 則可以得到星系模板庫的前m個(gè)主成分方向和超新星模板庫的前n個(gè)主成分方向分別為{gu},{sv},u= 1,2,··· ,m,v= 1,2,··· ,n.本文將主成分方向稱為光譜特征.

其中, 由于Nugent先生提供的模板在不同類別的超新星之間, 光譜數(shù)量存在較大差異. 由于其中一類最多只有12條光譜數(shù)據(jù), 于是對每一類在時(shí)間上等距選取12條光譜數(shù)據(jù)以用于構(gòu)建超新星的PCA模板庫.

圖3自上而下展示了星系與超新星模板庫按方差貢獻(xiàn)率排序后的前8個(gè)光譜特征. 此時(shí)選取超新星光譜特征累計(jì)貢獻(xiàn)率為95%的前8個(gè)光譜特征, 而對于星系成分選取累計(jì)貢獻(xiàn)率為90%的前180個(gè)光譜特征. 此時(shí), 可以將實(shí)測光譜看作是各光譜特征線性組合的近似[10,18], 即

圖3 光譜特征. 自上向下為星系(左)與超新星(右)模板庫按方差貢獻(xiàn)率排序后的前8個(gè)光譜特征, 為了方便展示, 對每條光譜特征進(jìn)行了偏移, 圖例數(shù)字表示偏移量.Fig.3 Spectral feature. The galaxy (left) and supernova (right) templates are sorted by variance contribution and the first eight spectral feature are top down displayed in figure. In order to display conveniently, each spectral feature has an offset, and the numbers in legend represent the value of offset.

設(shè)第k個(gè)混合光譜的系數(shù)向量ai=(a1,i,a2,i,··· ,am+n,i)T, 記

由(1)式有A·Ac=MT. 因此求解Ac即可求得各系數(shù)向量, 并實(shí)現(xiàn)光譜分解, 即計(jì)算

當(dāng)A不滿足正交性時(shí), 則可以用如下方法計(jì)算得到[10]:

其中, (β1,β2,··· ,βm,βm+1,··· ,βm+n)為A經(jīng)施密特正交化后所得的單位正交向量組.

通過光譜分解, 我們獲得了將一條光譜數(shù)據(jù)映射到以ai為坐標(biāo)向量的新特征空間的方法, 圖4分別展示了星系與超新星的分解效果. 圖5展示了部分變量組合的3維分布圖.

4.2 特征分布的相似性度量

假設(shè)已知兩條連續(xù)的概率密度曲線p1(x)、p2(x), 對p1、p2以?x為步長等距離離散可分別得到向量當(dāng)?x →0時(shí), 他們能夠很好地表征出該曲線的特征. 故定義兩條概率密度曲線的相似度:

圖4 星系光譜(左)與超新星光譜(右)分解效果, 圖中藍(lán)色曲線為經(jīng)過數(shù)據(jù)預(yù)處理的光譜, 紅色曲線是分解出來的超新星成分, 橙色曲線是分解出來的星系成分.Fig.4 Decomposition of Gal spectrum (left) and SN spectrum (right). The blue curve is the spectrum after data preprocessing, the red curve is the supernova component, and the orange curve is the galaxy component of the decomposition.

圖5 部分變量組合的3維分布圖Fig.5 Three dimensional distributions of partial variables combination

圖6為部分特征系數(shù)的頻率分布直方圖, 從圖中可以看到部分特征系數(shù)出現(xiàn)冗余的現(xiàn)象, 這些特征系數(shù)的分布在兩個(gè)類別之間的區(qū)別并不明顯, 甚至接近重合. 這些特征系數(shù)的存在對分類任務(wù)并沒有益處甚至?xí)Ψ诸愋Ч斐梢欢ǖ母蓴_.

于是需要度量一個(gè)特征系數(shù)的分布在兩個(gè)類別之間差異的不明顯程度, 即特征系數(shù)的分布相似程度, 以進(jìn)一步篩選特征系數(shù). 而(3)式恰好可以滿足這一需求.

圖7為相似度分布直方圖. 于是, 只要設(shè)定閾值s0, 保留Sims0的特征系數(shù), 即可過濾低效特征.

圖6 部分特征系數(shù)的頻率分布直方圖Fig.6 Frequency distribution histogram of partial feature coefficients

圖7 相似度分布直方圖Fig.7 Histogram of similarity

5 光譜識別

由圖5可見, 在光譜分解得到的特征空間下, 光譜識別問題可近似為一個(gè)線性可分的分類模型. 對此, 本文采用了線性支持向量機(jī)進(jìn)行分類訓(xùn)練.

5.1 支持向量機(jī)

線性可分是指在特征空間中能夠找到一個(gè)由下式表示的超平面來正確地描述分類任務(wù):

其中,w為超平面的法向量,x為坐標(biāo)向量. 在光譜分類問題中, 用第i條預(yù)處理光譜Mi提取的特征ai來代表該待分類光譜, 即有x=ai.

于是一個(gè)樣本點(diǎn)x到超平面的距離為:

若該超平面能夠正確分類, 即對于任一樣本(xi,yi), yi ∈{?1,1}, 當(dāng)xi為星系光譜時(shí),當(dāng)xi為超新星光譜時(shí)

只需要一定的事先變換, 樣本集總可以滿足:

此時(shí), 根據(jù)(5)式, 由該超平面所確定的分類間隔為2/||w||. 對于分類問題, 我們希望這個(gè)分類間距盡可能大, 而這等效于最小化因此我們只需要求解:

引入拉格朗日乘子ηi, 由拉格朗日乘子法有對偶問題

考慮防止過擬合或不是完全線性可分的情況下, 引入“松弛變量”ξi0,求解如下規(guī)劃問題:

即“軟間隔支持向量機(jī)”[12], 其中,C為常數(shù), 表示約束程度,C越大時(shí)松弛變量要越小, 通常取C=1.

在上文的特征提取方法下, 選取數(shù)據(jù)集的75%作為訓(xùn)練集, 剩余作為測試集. 迭代選取不同s0進(jìn)行SVM的訓(xùn)練, 得到了如圖8的誤判率變化曲線, 并在表3中展示了3個(gè)模型A、B、C的詳細(xì)效果.

表3 SVM在測試集上的表現(xiàn)Table 3 The SVM performance on the test set

圖8 誤判率曲線, s0為相似度閾值Fig.8 Mistake rate curve, s0 is the threshold of the similarity Sim

顯然, 相似度過濾算法能夠進(jìn)一步提高降維結(jié)果的質(zhì)量, 當(dāng)s0選取適當(dāng)時(shí), 仍能保證分類效果而數(shù)據(jù)維數(shù)較少. 在實(shí)際情況中, 可以視條件選擇合適的模型. 比如, 當(dāng)設(shè)備計(jì)算資源有限時(shí)采用算法空間復(fù)雜度較低的A模型, 當(dāng)計(jì)算資源足夠且希望漏判率盡可能小則使用漏判率最小的B模型.

由(5)式可以得到對于樣本x到?jīng)Q策面的距離r(x). 圖9展示了采用B模型時(shí)所得距離分布密度圖. SVM可以視為一個(gè)以原點(diǎn)為決策點(diǎn)的分類模型, 在直線L1右側(cè)即判斷為正, 在L1左側(cè)判斷為負(fù).

圖9 采用B模型時(shí), 樣本點(diǎn)到SVM超平面的距離r的分布. f1為星系樣本集上的距離分布曲線, f2為超新星樣本集上的距離分布曲線, L1和L2分別為SVM、NPSVM的決策線.Fig.9 Distribution of distance r from the sample point to the SVM hyperplane when using the model B.f1 is the distance distribution on the set of galaxies, f2 is the distance distribution on the set of supernovae. L1,L2 are the decision-making lines of SVM and NPSVM, respectively.

當(dāng)樣本數(shù)足夠大時(shí), 決策距離r的頻率分布曲線能夠近似于概率密度曲線, 即f1→p1(r)、f2→p2(r), 其中,f1為星系樣本集上r的頻率分布曲線,f2為超新星樣本集上r的頻率分布曲線;p1(r)為當(dāng)樣本x來自星系總體時(shí),r=r(x)在實(shí)數(shù)域R上的概率密度曲線,p2(r)為x來自超新星總體時(shí)r的概率密度曲線.

此時(shí), 設(shè)P1[r(x)]是樣本x來自星系總體的概率,P2[r(x)]為樣本x來自超新星總體的概率. 概率P1(r),P2(r)分別定義如下:

且給定x時(shí), 有P1(r)+P2(r)=1.

因此當(dāng)樣本量足夠大時(shí), 用f1(r)近似p1(r), 用f2(r)近似p2(r), 即可通過SVM模型計(jì)算一個(gè)樣本屬于兩個(gè)類別的概率. 當(dāng)有多個(gè)超新星候選體待觀測時(shí), 可以將候選體按照超新星的分類概率P1(r)排序, 制定觀測計(jì)劃, 從而更高效地尋找超新星.

5.2 NPSVM分類模型

來自超新星樣本集的樣本如果被誤判為星系, 即為超新星的漏判. 對于超新星巡天, 有時(shí)希望超新星的漏判率盡可能低. 為了進(jìn)一步降低漏判率, 不妨將距離r看作是經(jīng)過SVM映射所得的新特征.

設(shè)二元?jiǎng)澐諶= (R1,R2), 當(dāng)對于一個(gè)新的樣本(xi,yi),r(xi)經(jīng)過劃分R后判斷為星系時(shí), 記作r ∈R1; 當(dāng)r(xi)經(jīng)過R后判斷為超新星時(shí), 記為r ∈R2.

當(dāng)xi來自星系總體時(shí), 經(jīng)過劃分R后判斷為超新星則是誤判, 假設(shè)此事件發(fā)生的概率為Pe1; 同理, 對于超新星總體, 有Pe2. 由概率論有:

為了進(jìn)一步降低超新星漏判, 對于期望的超新星漏判上限ε0, 希望有Pe2=ε0, 此時(shí)要尋找滿足如下條件的最優(yōu)劃分R, 即Neyman-Pearson決策:

為了求解該問題, 可以構(gòu)建拉格朗日函數(shù)γ=Pe1+η(Pe2?ε0), 又因?yàn)??則有

其中,η為引入的拉格朗日乘子.

由拉格朗日乘子法有:

其中,r0由解出.

于是, 由(11)式和(13)式有最優(yōu)劃分:

其中,p1(r)/p2(r)稱為似然比. 于是, 通過調(diào)整ε0即可適應(yīng)不同的需求.

當(dāng)樣本數(shù)N足夠大時(shí), 有f1→p1(r)、f2→p2(r). 另外, 在實(shí)際計(jì)算中, 由于采樣的不足, 可能出現(xiàn)f1(r) = 0和f2(r) = 0的情況, 為了避免分母為0,頻率近似后可對決策函數(shù)(14)式做變換, 于是有:

其中,h(r)=f1(r0)f2(r)?f1(r)f2(r0),r0滿足本文稱(15)式為NPSVM模型.

若令ε0=0.002, 則得到圖9中的直線L2, 并計(jì)算h(r)得到其頻數(shù)分布圖如圖10,此時(shí)與SVM模型的比較效果如表4.

圖10 h(r)的分布直方圖Fig.10 Histogram of h(r)

表4 SVM與NPSVM效果比較Table 4 Comparison between SVM and NPSVM

計(jì)算得到超新星此時(shí)的漏判率為0.116%. 顯然, 效果比預(yù)期的ε0= 0.2%要好很多.但是, 星系被判為超新星的頻率為0.920%, 而該值采用普通的SVM時(shí), 僅有0.356%; 整體誤判率也由0.498%上升到0.854%. 故而, 此方法應(yīng)當(dāng)根據(jù)實(shí)際情況中的效益和需要適當(dāng)調(diào)整期望的漏判上限ε0.

6 總結(jié)與展望

本文通過數(shù)據(jù)預(yù)處理、特征提取、分類模型建立和優(yōu)化, 發(fā)現(xiàn)了光譜分解方法的冗余特征并給出了基于余弦相似度的優(yōu)化方案, 得到了一個(gè)效果良好的超新星識別模型. 相對于其他研究, 僅采用SVM模型時(shí), 本文能夠由決策距離分布給出一個(gè)樣本的分類概率, 以方便設(shè)置超新星的后續(xù)觀測優(yōu)先級. 而追求較小超新星漏判率時(shí), 本文提出的NPSVM模型也能夠有效限制其漏判率.

但是在實(shí)踐中, 雖然NPSVM模型可以降低漏判率, 但也會增加將星系判為超新星的概率. 超新星的最終證認(rèn)往往還要通過人工證認(rèn)或后續(xù)觀測進(jìn)行確認(rèn), 而星系被判為超新星會增加人力與時(shí)間成本, 因而需要根據(jù)實(shí)際情況進(jìn)行漏判率的限制. 如果希望整體誤判率最小, 還可以選擇僅采用SVM模型而不采用Neyman-Pearson決策.

值得一提, 本文所涉及的PCA、SVM算法, 當(dāng)數(shù)據(jù)量越大時(shí), 在統(tǒng)計(jì)意義下越可靠,效果也越好. 預(yù)期隨著更多的超新星發(fā)現(xiàn), 模型效果將不斷提高, 從而形成正反饋機(jī)制推動超新星的發(fā)現(xiàn).

致謝感謝審稿人對文章提出的寶貴建議, 使得文章的質(zhì)量有了顯著的提高. 感謝WISeREP、CfA、SDSS公開的光譜數(shù)據(jù)和Nugent公開的超新星模板和相關(guān)工作者的付出.