李 慧 丁雨君 李鄉(xiāng)儒 張金區(qū)

(1 華南師范大學(xué)計(jì)算機(jī)學(xué)院 廣州 510631)

(2 華南師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院 廣州 510631)

1 引言

射電天文學(xué)是現(xiàn)代天文學(xué)的重要分支[1–2]. 由于儀器的高靈敏度和電子類(lèi)產(chǎn)品的普及, 射電天文數(shù)據(jù)會(huì)受到人類(lèi)生產(chǎn)、生活以及非觀測(cè)目標(biāo)的宇宙輻射源產(chǎn)生的射電信號(hào)的影響, 此類(lèi)影響稱(chēng)之為射頻干擾(Radio Frequency Interference, RFI)[3–5].宇宙輻射源, 如太陽(yáng), 可能會(huì)干擾觀測(cè). 這類(lèi)干擾由于位置和強(qiáng)度已知, 易于回避或糾正. 但人為活動(dòng)產(chǎn)生的干擾, 通常不可預(yù)測(cè)且不穩(wěn)定、強(qiáng)度未知、難以控制[6]. 常見(jiàn)的射頻干擾源包括電視信號(hào)、調(diào)頻無(wú)線電傳輸、全球定位系統(tǒng)、手機(jī)和飛機(jī)導(dǎo)航通訊等[7]. 不同射頻干擾源的頻率和時(shí)間特性有差異, 導(dǎo)致整個(gè)RFI檢測(cè)問(wèn)題很復(fù)雜. 射頻干擾檢測(cè)成為射電觀測(cè)數(shù)據(jù)處理中遇到的重要挑戰(zhàn)之一.

因此, 該問(wèn)題備受關(guān)注, 而且研究者給出了一系列RFI檢測(cè)方法. 從原理上來(lái)說(shuō), RFI檢測(cè)方法大致可分為成分分解法、閾值分析法和機(jī)器學(xué)習(xí)法等. 成分分解法的基本思想是從數(shù)據(jù)中自動(dòng)發(fā)現(xiàn)RFI在時(shí)間或頻率方面展現(xiàn)出來(lái)的規(guī)律性, 據(jù)此實(shí)現(xiàn)對(duì)RFI與非RFI數(shù)據(jù)成分的分離. 這類(lèi)方法適用于RFI在時(shí)間或頻率上表現(xiàn)出重復(fù)模式的情況, 但是不能處理各種不規(guī)則信號(hào)[1]. 例如, 奇異值分解法(Singular Value Decomposition,SVD)[8]和主成分分析法(Principle Component Analysis, PCA)[9–10]都是典型的成分分解類(lèi)RFI檢測(cè)方法.

隨著機(jī)器學(xué)習(xí)的迅猛發(fā)展和廣泛應(yīng)用, 聚類(lèi)分析法、卷積神經(jīng)網(wǎng)絡(luò)法等機(jī)器學(xué)習(xí)方法在RFI檢測(cè)[11–13]中的應(yīng)用漸受關(guān)注. 機(jī)器學(xué)習(xí)算法可從海量數(shù)據(jù)中學(xué)習(xí)數(shù)據(jù)的知識(shí)表示和數(shù)據(jù)成分之間的復(fù)雜關(guān)系,進(jìn)而完成模式的發(fā)現(xiàn)或識(shí)別,并在RFI檢測(cè)中初步展示出了不錯(cuò)的應(yīng)用潛力. 但是機(jī)器學(xué)習(xí)模型訓(xùn)練耗時(shí), 模型正確性驗(yàn)證復(fù)雜, 而且樣本數(shù)據(jù)和特征的選擇直接影響其分類(lèi)精度, 特別是有標(biāo)注樣本數(shù)據(jù)不足時(shí)難以得到高精度的訓(xùn)練模型.

閾值分析法因其實(shí)現(xiàn)簡(jiǎn)單、檢測(cè)結(jié)果精度較高而被廣泛應(yīng)用. 它的理論依據(jù)是射電天文望遠(yuǎn)鏡所接收到的來(lái)自地球輻射源產(chǎn)生的RFI信號(hào)的強(qiáng)度往往大于來(lái)自太空的信號(hào)強(qiáng)度. 因此, 當(dāng)一個(gè)信號(hào)強(qiáng)度值超過(guò)某個(gè)閾值時(shí), 閾值分析法可將它標(biāo)記為RFI. 代表性閾值分析法有CUSUM (cumulative sum)法[14]、Simple Thresholding法[15]和Combinatorial Thresholding法[16]等. CUSUM算法通過(guò)估算累積樣本的方差或平均值得到閾值, 高于閾值的數(shù)據(jù)被標(biāo)記為RFI. 此方法簡(jiǎn)單、快速, 但是無(wú)法檢測(cè)到RFI出現(xiàn)的準(zhǔn)確起始時(shí)間, 只能用于估計(jì)其粗略范圍. 因此, CUSUM適合于RFI的預(yù)檢, 即首先用CUSUM發(fā)現(xiàn)存在RFI的粗略位置, 然后用其他更準(zhǔn)確但相對(duì)耗時(shí)的方法精確標(biāo)記[16]. Simple Thresholding方法在檢測(cè)某行(列)的觀測(cè)數(shù)據(jù)時(shí), 使用該行(列)的中位數(shù)作為閾值.該方法運(yùn)行速度快, 但在檢測(cè)瞬時(shí)射頻干擾時(shí),僅可以檢測(cè)到干擾的峰值, 通常會(huì)忽略上升過(guò)程中部分強(qiáng)度弱的RFI樣本點(diǎn). 然而, RFI干擾往往會(huì)影響相鄰位置的多個(gè)樣本. 因此, Simple Thresholding方法易造成漏檢. 為此, Combinatorial Thresholding算法通過(guò)采用滑動(dòng)窗口與多次迭代的機(jī)制進(jìn)行檢測(cè), 當(dāng)樣本組合的值超過(guò)閾值時(shí), 將該樣本組合標(biāo)記為RFI. Combinatorial Thresholding算法解決了Simple Thresholding算法瞬時(shí)射頻干擾RFI漏檢的問(wèn)題, 但存在RFI過(guò)檢測(cè)傾向[15]. 因此, 學(xué)者們提出了改進(jìn)算法VarThreshold和SumThreshold. 當(dāng)檢測(cè)窗口內(nèi)所有樣本的讀數(shù)均大于閾值時(shí), VarThreshold算法將窗口內(nèi)的樣本標(biāo)記為RFI. 而SumThreshold算法則僅關(guān)注窗口內(nèi)尚未被檢測(cè)為RFI的像素: 如果這些像素的均值大于給定閾值, 則將它們標(biāo)記為RFI. RFI強(qiáng)度可變、形態(tài)多樣以及不可預(yù)測(cè)的性質(zhì)使得檢測(cè)具有挑戰(zhàn)性, 構(gòu)建穩(wěn)健的RFI檢測(cè)方法至關(guān)重要[17]. 研究表明, SumThreshold方法在RFI檢測(cè)中具有較高的精度[16], 已經(jīng)廣泛地應(yīng)用于射電天文數(shù)據(jù)處理中[18–21], 成為RFI檢測(cè)的典型算法. 因此, 本文從原理、性能、優(yōu)化等方面對(duì)SumThreshold算法進(jìn)行深入探討, 以期促進(jìn)Sum-Threshold算法的研究和應(yīng)用.

2 SumThreshold算法描述

在射電天文學(xué)中, 射頻干擾可分為3類(lèi): 脈沖干擾、長(zhǎng)窄帶干擾和復(fù)合干擾[15]. 脈沖干擾是指在短時(shí)間內(nèi)出現(xiàn)的較寬頻帶干擾. 長(zhǎng)窄帶干擾中干擾是出現(xiàn)在某個(gè)小的頻率子帶內(nèi)的一個(gè)相對(duì)恒定的流. 復(fù)合干擾是前兩種干擾的組合. 由此可知, RFI干擾在時(shí)間或頻率上會(huì)影響多個(gè)位置連續(xù)的像素.傳統(tǒng)閾值分析法是通過(guò)對(duì)單一數(shù)據(jù)值與某個(gè)指定的閾值做比較實(shí)現(xiàn)RFI檢測(cè). 但是, 這類(lèi)基于單個(gè)像素比較的閾值類(lèi)方法也存在一定的局限性: 會(huì)引起個(gè)別像素的RFI假陽(yáng)性或假陰性. 對(duì)此, 一種解決思路是將像素鄰近性因素納入考慮之中.

為了充分利用RFI數(shù)據(jù)間的位置鄰近性, Sum-Threshold算法構(gòu)造了滑動(dòng)窗口和閾值集合. 當(dāng)檢測(cè)窗口內(nèi)數(shù)據(jù)的均值高于閾值時(shí), 窗口內(nèi)數(shù)據(jù)被標(biāo)記為RFI. 通過(guò)窗口的滑動(dòng)和迭代, 實(shí)現(xiàn)對(duì)整個(gè)觀測(cè)數(shù)據(jù)的RFI檢測(cè). SumThreshold的算法流程如圖1所示.

圖1 SumThreshold的算法流程圖Fig.1 The flow chart of SumThreshold method

2.1 基線矯正

閾值法進(jìn)行RFI檢測(cè)的基本前提是: 如果數(shù)據(jù)沒(méi)有受到RFI干擾且未疊加射電天文信號(hào), 則相應(yīng)的讀數(shù)基本恒定[22], 并將其簡(jiǎn)稱(chēng)為背景響應(yīng)恒定假設(shè). 理想情況下, 基線在頻域和時(shí)域中保持恒定. 但是, 幾乎所有的射電天文數(shù)據(jù)受到系統(tǒng)漂移、大氣效應(yīng)、地面輻射等因素的影響, 由此導(dǎo)致背景恒定假設(shè)一般不成立. 這種背景不恒定性在有些文獻(xiàn)中稱(chēng)為基線變化. 例如, 圖2展示了一幅500 m口徑球面射電望遠(yuǎn)鏡(Five-hundred-meter Aperture Spherical radio Telescope,FAST)的時(shí)間-頻率觀測(cè)圖像[23], 觀測(cè)時(shí)間為0.05 s, 頻率范圍是1000–1500 MHz; 在該觀測(cè)中基線變化導(dǎo)致低頻帶像素讀數(shù)較低, 且像素讀數(shù)隨著頻率的增大而整體上增大. 這種基線變化致使非RFI數(shù)據(jù)讀數(shù)取值區(qū)間范圍增大, 無(wú)法直接通過(guò)某個(gè)固定閾值實(shí)現(xiàn)RFI檢測(cè).為此, 在運(yùn)用SumThreshold算法檢測(cè)RFI之前, 需進(jìn)行基線擬合和剔除, 盡量減少背景不恒定性造成的影響.

圖2 基線對(duì)閾值類(lèi)方法的影響. 這是FAST的一個(gè)時(shí)間-頻率觀測(cè)圖像, 觀測(cè)時(shí)間為0.05 s, 頻率范圍是1000–1500 MHz[23]. (a)原始觀測(cè)數(shù)據(jù);(b)去除基線效應(yīng)后的結(jié)果. 原始觀測(cè)數(shù)據(jù)中不同頻帶像素讀數(shù)浮動(dòng)范圍變化非常大, 對(duì)RFI檢測(cè)中的閾值設(shè)定造成困擾.Fig.2 The influence of baseline on the threshold method. This is a time-frequency image observed using FAST for 0.05 s in the frequency range from 1000 to 1500 MHz[23]. (a) An original time-frequency image; (b) the result after baseline removal. The range of pixel intensity is very broad in different frequency bands of the original observation, so that it is hard to set an appropriate threshold for RFI detection.

文獻(xiàn)[23]運(yùn)用非對(duì)稱(chēng)加權(quán)懲罰最小二乘法(Asymmetrically reweighted Penalized Least Squares, ArPLS)進(jìn)行基線擬合和剔除. 與傳統(tǒng)的二維低階多項(xiàng)式法[16]相比, 此方法更加高效、準(zhǔn)確和穩(wěn)健. 基線剔除修正了原圖數(shù)據(jù)的背景不一致性, 原始圖像中對(duì)比度低的區(qū)域變得易于檢測(cè)RFI (圖2).

2.2 RFI檢測(cè)方向指定

對(duì)于時(shí)間-頻率觀測(cè)數(shù)據(jù), 基于SumThreshold算法的RFI檢測(cè)有3個(gè)實(shí)施策略: 基于時(shí)間維度的RFI檢測(cè)、基于頻率維度的RFI檢測(cè)以及基于時(shí)間和頻率的雙向檢測(cè). 基于時(shí)間和頻率的雙向檢測(cè), 首先沿著時(shí)間維度檢測(cè)RFI, 然后基于第1次檢測(cè)結(jié)果,在頻率方向繼續(xù)進(jìn)行RFI檢測(cè),這種檢測(cè)策略簡(jiǎn)稱(chēng)為時(shí)間-頻率雙向檢測(cè), 反之則稱(chēng)為頻率-時(shí)間雙向檢測(cè). 為了比較3種實(shí)施策略的效果差異, 本文基于FAST觀測(cè)數(shù)據(jù)[23]進(jìn)行實(shí)驗(yàn). 圖3的FAST觀測(cè)數(shù)據(jù)中由于存在取值較大的RFI, 導(dǎo)致圖像細(xì)節(jié)不明顯. SumThreshold算法單一方向RFI檢測(cè)結(jié)果對(duì)比見(jiàn)圖4. 圖4 (a)、(b)為沿時(shí)間維度或頻率維度的單一方向RFI檢測(cè)結(jié)果. 圖4 (c)、(d)展示的是僅一方向檢測(cè)發(fā)現(xiàn)的RFI. 實(shí)驗(yàn)結(jié)果表明, 兩種單一方向的檢測(cè)策略實(shí)驗(yàn)結(jié)果存在差異. 但是, 強(qiáng)干擾在兩種檢測(cè)策略里均被正確檢測(cè), 而存在差異的是較弱的RFI.

圖3 FAST觀測(cè)數(shù)據(jù). 觀測(cè)時(shí)間為0.05 s, 頻率范圍是1000–1500 MHz. 由于存在取值較大的RFI, 導(dǎo)致圖像細(xì)節(jié)不明顯.該數(shù)據(jù)來(lái)源于FAST觀測(cè)[23].Fig.3 A time-frequency image of FAST observation for 0.05 s in the frequency range from 1000 to 1500 MHz. There are some RFI with extremely large values. Therefore,visibility of the image details is poor. The experimental data are of a FAST observation[23].

因此, 在實(shí)際應(yīng)用中一般采用基于時(shí)間和頻率的雙向檢測(cè)策略. 圖5 (a)是先沿著時(shí)間方向再沿著頻率方向的RFI檢測(cè)結(jié)果, 圖5 (b)是先沿著頻率方向再沿著時(shí)間方向的RFI檢測(cè)結(jié)果. 與單一方向檢測(cè)結(jié)果(圖4)相比, 雙向RFI檢測(cè)能夠更全面地發(fā)現(xiàn)RFI. 圖6為基于不同雙向檢測(cè)策略的SumThreshold算法效果對(duì)比, 結(jié)果表明頻率和時(shí)間檢測(cè)方向的先后順序不同時(shí), RFI檢測(cè)結(jié)果有一定差異.基于雙向檢測(cè)和單向檢測(cè)的SumThreshold算法檢測(cè)結(jié)果比較如圖7所示.將單一時(shí)間維度RFI檢測(cè)結(jié)果和頻率維度RFI檢測(cè)結(jié)果合并, 與雙向RFI檢測(cè)結(jié)果對(duì)比發(fā)現(xiàn), 前者檢測(cè)RFI數(shù)量大于后者. 這主要是因?yàn)? 基于雙向標(biāo)記策略的SumThreshold算法在標(biāo)記過(guò)程中, 當(dāng)?shù)?個(gè)方向完成檢測(cè)時(shí), 會(huì)對(duì)下一個(gè)方向的檢測(cè)產(chǎn)生影響, 即部分樣本點(diǎn)在上一個(gè)方向的檢測(cè)中已經(jīng)被標(biāo)記, 從而導(dǎo)致在進(jìn)行下一方向檢測(cè)時(shí), 未標(biāo)記樣本點(diǎn)的平均值減小. 在閾值集合不變的情況下, 第2方向檢測(cè)的RFI數(shù)量減少. 因此, 可以在第2方向檢測(cè)時(shí), 根據(jù)數(shù)據(jù)的分布特征,調(diào)整閾值參數(shù)χ1, 進(jìn)而得到更優(yōu)的閾值集合(閾值集合的詳細(xì)討論, 見(jiàn)本文2.3節(jié)).

圖4 SumThreshold (ST)算法的不同實(shí)施策略效果比較. (a) ST算法沿時(shí)間維度的RFI檢測(cè)結(jié)果; (b) ST算法沿頻率維度的RFI檢測(cè)結(jié)果;(c)在圖(a)中被檢測(cè)到但未在圖(b)中檢測(cè)到的RFI; (d)在圖(b)中被檢測(cè)到但未在圖(a)中檢測(cè)到的RFI. 圖中白色點(diǎn)表示RFI數(shù)據(jù). 由于強(qiáng)干擾被兩種策略均正確檢測(cè), 所以在圖(c)和(d)中均未標(biāo)記, 且它們比圖(a)和圖(b)整體上顯得暗一些. 該數(shù)據(jù)來(lái)源于FAST觀測(cè)[23].Fig.4 The detection results of SumThreshold with different implementation schemes. (a) The results detected along time direction; (b) the results detected along frequency direction; (c) the RFI detected only along time direction, but not frequency direction; (d) the RFI detected only along frequency direction, but not time direction. The RFI data are presented in white.Strong RFI are successfully detected by both strategies. Therefore, these strong RFI are not showed in panels (c) and (d), and panels (c) and (d) are darker than panels (a) and (b). The experimental data are of a FAST observation[23].

圖5 雙向SumThreshold算法的RFI檢測(cè)結(jié)果. (a)先沿著時(shí)間方向再沿著頻率方向的RFI檢測(cè)結(jié)果(簡(jiǎn)稱(chēng)為時(shí)間-頻率雙向檢測(cè)); (b)先沿著頻率方向再沿著時(shí)間方向的RFI檢測(cè)結(jié)果(簡(jiǎn)稱(chēng)為頻率-時(shí)間雙向檢測(cè)); (c)將圖4 (a)和(b)單一方向檢測(cè)結(jié)果合并后的RFI標(biāo)記結(jié)果. 該數(shù)據(jù)來(lái)源于FAST觀測(cè)[23].Fig.5 The RFI detection results based on bidirection detection. (a) The RFI detection results detected firstly along time direction and then frequency direction (time-frequency bidirection detection); (b) the RFI detection results firstly along frequency direction and then time direction (frequency-time bidirection detection); (c) the RFI detection result from panels 4(a) and (b). The experimental data are of a FAST observation[23].

圖6 基于不同雙向檢測(cè)策略的SumThreshold算法效果對(duì)比. (a)被時(shí)間-頻率雙向方式檢測(cè)到但未被頻率-時(shí)間雙向檢測(cè)方式發(fā)現(xiàn)的RFI; (b)被頻率-時(shí)間雙向檢測(cè)方式發(fā)現(xiàn)但未被時(shí)間-頻率雙向檢測(cè)方式發(fā)現(xiàn)的RFI. 由于強(qiáng)干擾被兩種雙向策略均正確檢測(cè), 所以它們?cè)谶@個(gè)實(shí)驗(yàn)中均未標(biāo)記, 且該圖比圖5整體上暗一些. 該數(shù)據(jù)來(lái)源于FAST觀測(cè)[23].Fig.6 Detection result comparison of the SumThreshold with different bidirection detection schemes. (a) The RFI detected by SumThreshold with time-frequency scheme; (b) the RFI detected by SumThreshold with frequency-time scheme. Strong RFI are successfully detected by both strategies. Therefore, the strong RFI are not showed in Fig.5 but exist in (a) and (b) of Fig.6.The existences of the strong RFI make the panels (a) and (b) of Fig.6 darker than Fig.5. The experimental data are of a FAST observation[23].

2.3 使用SumThreshold算法對(duì)RFI進(jìn)行檢測(cè)

使用SumThreshold算法對(duì)RFI做檢測(cè)需要進(jìn)行多次迭代. 每次迭代, 使用某個(gè)指定大小的窗口沿著時(shí)間或頻率方向進(jìn)行移動(dòng), 并對(duì)窗口內(nèi)的像素求平均值, 據(jù)此進(jìn)行閾值檢測(cè). 檢測(cè)窗口的尺寸隨著迭代次數(shù)的增加而增大. 例如, 當(dāng)時(shí)頻數(shù)據(jù)有N個(gè)觀測(cè)頻帶時(shí), 跨頻帶方向移動(dòng)檢測(cè)的情況下最多可迭代N次. 實(shí)際應(yīng)用中, 一般不需要進(jìn)行最大次數(shù)的迭代即可取得良好結(jié)果. 文獻(xiàn)[23]中根據(jù)FAST數(shù)據(jù)的特點(diǎn), 經(jīng)驗(yàn)性地將迭代次數(shù)設(shè)置為11.

2.3.1 RFI數(shù)據(jù)標(biāo)記準(zhǔn)則

在SumThreshold算法每次迭代中, 給定迭代窗口的大小和相應(yīng)的檢測(cè)閾值. 根據(jù)算法原理可知, 之前迭代中已被標(biāo)記RFI的像素在新一輪迭代中將被替換為當(dāng)前的檢測(cè)閾值. 這樣的數(shù)據(jù)重置,可避免一些RFI像素造成周?chē)荝FI像素被錯(cuò)誤檢測(cè)為RFI(假陽(yáng)性). 因此, 對(duì)于每個(gè)檢測(cè)窗口, 實(shí)際上是計(jì)算當(dāng)前窗口內(nèi)尚未被檢測(cè)為RFI像素的均值, 據(jù)此判斷這些像素是否為RFI成分[16,21,23]. 基于SumThreshold的RFI檢測(cè)原理可形式化表達(dá)為:

(1)式中第p+1次迭代, 滑動(dòng)窗口大小為正整數(shù)Mp+1, 檢測(cè)閾值是χp+1.Ri是下標(biāo)為i的像素?cái)?shù)據(jù),i是數(shù)據(jù)下標(biāo).表示第p+1次迭代中,以下標(biāo)為v的數(shù)據(jù)作為起點(diǎn), 窗口大小為Mp+1的滑動(dòng)窗口標(biāo)記集合, 設(shè)如果Rv被檢測(cè)為RFI, 則fpv= 1, 否則,fpv= 0.fpv取值為1表示p次迭代后此數(shù)據(jù)已標(biāo)記為RFI數(shù)據(jù), 反之, 此數(shù)據(jù)在當(dāng)前迭代步驟被判定為非RFI數(shù)據(jù).Count是窗口內(nèi)尚未被標(biāo)記為RFI像素的個(gè)數(shù):

閾值χp+1是在VarThreshold算法閾值公式基礎(chǔ)上,通過(guò)參數(shù)的優(yōu)化來(lái)獲取[16]. VarThreshod算法閾值公式為:

經(jīng)驗(yàn)表明,ρ= 1.5時(shí)SumThreshold算法具有較好的效果. 為了確定閾值χ1, 可在某個(gè)給定的觀測(cè)數(shù)據(jù)集上最小化錯(cuò)誤概率實(shí)現(xiàn), 在這過(guò)程中參數(shù)ρ保持不變. 確定χ1之后, 根據(jù)(3)式可計(jì)算得到各χp+1,p≥1. 文獻(xiàn)[23]以及SEEK[21]軟件中對(duì)于SumThreshold算法的閾值計(jì)算, 設(shè)置了數(shù)組η.由(3)式計(jì)算出的值除以對(duì)應(yīng)的數(shù)組η值, 構(gòu)成最終的閾值集合.

文獻(xiàn)[18]對(duì)于RFI數(shù)據(jù)的判斷, 未采用當(dāng)前未標(biāo)記為RFI像素的均值, 而是計(jì)算像素絕對(duì)值的均值. 這種RFI標(biāo)記準(zhǔn)則適用于射頻干擾振蕩嚴(yán)重,但像素平均值卻接近于0的射電數(shù)據(jù). 因此, 應(yīng)當(dāng)具體問(wèn)題具體分析, 在不同應(yīng)用背景中, 選取適宜的RFI標(biāo)記準(zhǔn)則. 本文后續(xù)內(nèi)容仍然采用(1)式的RFI標(biāo)記準(zhǔn)則.

2.3.2 SumThreshold算法設(shè)計(jì)

對(duì)于某個(gè)給定寬度為Mn的檢測(cè)窗口, Sum-Threshold檢測(cè)的過(guò)程請(qǐng)見(jiàn)算法1:

算法1:

輸入: 被檢測(cè)數(shù)據(jù)長(zhǎng)度為L(zhǎng), 數(shù)據(jù)集為{R0,R1, ··· , RL-1}, 算法迭代總次數(shù)為max, 窗口大小的集合{M1,M2,··· ,Mmax}, 閾值集合{χ1, χ2,··· , χmax}, 之前迭代生成的標(biāo)記集合{f0, f1, ··· ,fL-1}, 標(biāo)記集合元素初始化值為0

輸出: 標(biāo)記集合{f0,f1,··· ,fL-1}

(1)獲取當(dāng)前窗口大小Mn、閾值χn、初始化變量i←0、未檢測(cè)為RFI的數(shù)據(jù)之和Sum←0、Count←0, 本次迭代生成的標(biāo)記集合{t0, t1, ··· , tL-1}←{f0,f1,··· ,fL-1}

(2) Whilei/=Mndo

(3) Iffi== 0 then

(4) Sum←Sum +Ri

(5) Count←Count + 1

(6) End if

(7)i ←i+1

(8) End While

(9) Whilei/=Ldo

(10) If Sum/Count＞χnor Sum/Count＜-χnthen

(11) Forj ∈{i-Mn,··· ,i-1}do

(12)tj ←1

(13) End For

(14) End if

(15) Iffi== 0 then

(16) Sum←Sum +Ri

(17) Count←Count + 1

(18) End if

(19) Iffi-Mn== 0 then

(20) Sum←Sum- Ri-Mn

(21) Count←Count-1

(22) End if

(23)i ←i+1

(24) End While

(25){f0,f1,··· ,fL-1}←{t0,t1,··· ,tL-1}

步 驟(2)–(6), 計(jì) 算 從 第1個(gè) 數(shù) 據(jù) 開(kāi) 始, 長(zhǎng)為Mn的窗口內(nèi)尚未被檢測(cè)為RFI的像素的均值.步驟(10)–(14)判斷當(dāng)前窗口內(nèi)像素的均值是否大于閾值, 如果判斷結(jié)果為“是”, 則將當(dāng)前窗口內(nèi)所有尚未被檢測(cè)為射頻干擾的像素標(biāo)記為RFI. 步驟(15)–(18)的意思是, 如果檢測(cè)窗口外右側(cè)像素在之前的迭代中未被判斷為RFI數(shù)據(jù), 則將此數(shù)據(jù)移入窗口, 窗口內(nèi)數(shù)據(jù)之和以及數(shù)據(jù)個(gè)數(shù)增加. 步驟(19)–(22)的意思是如果當(dāng)前窗口內(nèi)左側(cè)像素在之前的迭代中未被判斷為RFI數(shù)據(jù), 則將其移出窗口, 窗口內(nèi)數(shù)據(jù)之和以及數(shù)據(jù)個(gè)數(shù)相應(yīng)減少. 循環(huán)執(zhí)行步驟(10)–(22), 直至所有數(shù)據(jù)被處理, 更新標(biāo)記集合, 本次迭代結(jié)束.

2.3.3 基于SumThreshold算法的RFI檢測(cè)示例

圖8顯示了一個(gè)觀測(cè)數(shù)據(jù)的截面及其RFI標(biāo)記結(jié)果的示例. 該數(shù)據(jù)是從某個(gè)頻帶中截取得到. 因此, 圖8 (a)的橫坐標(biāo)和縱坐標(biāo)分別表示時(shí)間和流量. 經(jīng)過(guò)基線擬合和剔除后, 各像素的值分別為:{1, 1, 2, 1, 3, 1, 12, 14, 16, 15, 17, 20, 24, 26, 22,18, 14, 13, 11, 1, 1, 3, 3, 1, 1, 1, 2, 3, 1, 1, 1, 1}.基于SumThreshold算法進(jìn)行RFI檢測(cè)后, 值為{12,14, 16, 15, 17, 20, 24, 26, 22, 18, 14, 13, 11}的像素被標(biāo)記為RFI (圖8 (b)中斜線標(biāo)記的部分).

圖8 一個(gè)觀測(cè)數(shù)據(jù)及其RFI檢測(cè)結(jié)果的截面. (a)一個(gè)觀測(cè)數(shù)據(jù)的截面; (b)左圖數(shù)據(jù)的RFI檢測(cè)結(jié)果.Fig.8 The intersecting surface of an observation and its RFI detection results. (a) An intersecting surface of an observation; (b)an intersecting surface of RFI detection results of the left observation.

圖8(a)的RFI檢測(cè)過(guò)程如圖9所示. 在該實(shí)驗(yàn)中, SumThreshold算法經(jīng)過(guò)了3次迭代, 第1次迭代滑動(dòng)窗口寬度是M1= 1, 閾值為χ1= 20. 迭代輸入數(shù)據(jù)為圖8 (a)中的數(shù)據(jù), 如圖9 (a)所示. 因?yàn)榇翱诖笮?, 隨著窗口滑動(dòng), 每一個(gè)數(shù)據(jù)均與閾值進(jìn)行比較, 值為24、26和22的數(shù)據(jù)大于閾值. 因此, 相應(yīng)像素被標(biāo)記為RFI (圖9 (b)). 第2次迭代的滑動(dòng)窗口寬度為M2= 2, 閾值是χ2= 13.3333.為了防止已標(biāo)記為RFI的數(shù)據(jù)值過(guò)大導(dǎo)致周?chē)南袼乇诲e(cuò)誤檢測(cè)為射頻干擾, 在第2次迭代檢測(cè)之前將第1輪迭代時(shí)被標(biāo)記為RFI的像素?cái)?shù)值改為當(dāng)前閾值13.3333 (圖9 (c)). 第2次迭代后新被標(biāo)記為RFI的數(shù)據(jù)分別為14、16、15、17、20、18、14和13 (圖9 (d)). 對(duì)于值為13的數(shù)據(jù), 當(dāng)其進(jìn)入滑動(dòng)窗口時(shí), 窗口內(nèi)的數(shù)據(jù)值為14和13, 該窗口內(nèi)像素的均值大于閾值13.3333, 因此標(biāo)記數(shù)據(jù)13為RFI. 如果將值為13的像素單獨(dú)與閾值相比,則由于小于閾值而不會(huì)被檢測(cè)到, 但是因?yàn)榇翱趦?nèi)像素的整體平均效應(yīng), 該像素被SumThreshold算法標(biāo)記為RFI, 解決了傳統(tǒng)閾值法RFI由單個(gè)像素比較導(dǎo)致的漏檢問(wèn)題. 第3次迭代中窗口寬度是M3= 4, 檢測(cè)閾值為χ3= 10.5180. 之前迭代中被標(biāo)記為RFI的像素的值在本次迭代中置換為10.5180 (圖9 (e)). 通過(guò)與第2次迭代相同操作流程, 滑動(dòng)窗口對(duì)數(shù)據(jù)處理后, 新被標(biāo)記為RFI的數(shù)據(jù)是值為12和11的像素(圖9 (f)).

圖9 RFI檢測(cè)過(guò)程示例圖. 圖(a)、(c)、(e)分別是第1、2、3次迭代前的輸入數(shù)據(jù); 圖(b)、(d)、(f)分別是第1, 2, 3次迭代處理后的數(shù)據(jù). 詳述見(jiàn)2.3節(jié).Fig.9 The iterative process of RFI detection. Panels (a), (c), (e) are the input data of the iterative process and panels (b), (d),(f) are the results of the iterative process. More can be found in Sec. 2.3.

2.4 算法性能分析

2.4.1 時(shí)間復(fù)雜度

假設(shè)待檢測(cè)數(shù)據(jù)長(zhǎng)度為L(zhǎng), 則理論上檢測(cè)窗口大小的集合為{1,2,··· ,L-1}. 每次迭代要求對(duì)每個(gè)窗口子序列大小的數(shù)據(jù)進(jìn)行傳遞. SumThreshold算法的時(shí)間復(fù)雜度為O(L2). 這對(duì)于算法效率要求較高的應(yīng)用很難滿足要求. 因此, 基于指數(shù)增長(zhǎng)設(shè)定窗口大小, 則檢測(cè)窗口集合成為原始窗口集合的子集. 例如, 若將窗口大小設(shè)為[1,2,4,8,16,···],則算法時(shí)間復(fù)雜度為O(Llog2L). 與理論上的檢測(cè)窗口集合相比, 減少了尺寸為[3,5,6,7,···]的窗口.這些被減少的窗口所檢測(cè)的數(shù)據(jù)特征, 仍可能被約簡(jiǎn)窗口集合檢測(cè)到. 因此, 上述對(duì)候選窗口減少的做法對(duì)算法精度不產(chǎn)生顯著影響.

在減少候選檢測(cè)窗口集合的基礎(chǔ)上, 設(shè)定數(shù)據(jù)集子集大小, 可進(jìn)一步將算法復(fù)雜度降低. 在文獻(xiàn)[24]中, 通過(guò)限定數(shù)據(jù)子集小于1024, 算法效率明顯提高. 但是算法的精度會(huì)受影響, 這是因?yàn)榉浅Ｎ⑷趸蛘叻秶蟮奶卣鲗⒈缓雎?

2.4.2 算法適用性

已有研究表明, SumThreshold算法是射電數(shù)據(jù)RFI檢測(cè)的有效方法, 適合所有類(lèi)型的RFI檢測(cè).根據(jù)各類(lèi)RFI的特性可知, 基于頻率方向的檢測(cè)適用于帶狀射頻干擾, 基于時(shí)間方向的檢測(cè)則對(duì)突發(fā)脈沖射頻干擾效果較好[15]. 基于時(shí)間和頻率雙向的檢測(cè), 可檢測(cè)復(fù)合干擾. 與其他閾值類(lèi)RFI檢測(cè)算法一樣, SumThreshold算法對(duì)帶狀RFI和復(fù)合型RFI檢測(cè)均受射電數(shù)據(jù)頻率子帶大小的影響[15].

2.4.3 算法精度

SumThreshold算法對(duì)射電數(shù)據(jù)進(jìn)行RFI標(biāo)記,結(jié)果可分為4類(lèi): 真陽(yáng)性(True positive, TP): 標(biāo)記為RFI的數(shù)據(jù), 實(shí)際上也是RFI數(shù)據(jù); 假陽(yáng)性(False positive, FP): 標(biāo)記為RFI的數(shù)據(jù), 實(shí)際上不是RFI數(shù)據(jù); 真陰性(True negative, TN): 標(biāo)記為非RFI的數(shù)據(jù), 實(shí)際上也是非RFI數(shù)據(jù); 假陰性(False negative, FN): 標(biāo)記為非RFI的數(shù)據(jù), 實(shí)際上是RFI數(shù)據(jù). 這4類(lèi)結(jié)果可構(gòu)成混淆矩陣, 見(jiàn)表1.

基于表1, 射頻干擾檢測(cè)結(jié)果的常用精度度量指標(biāo)有:

表1 RFI檢測(cè)混淆矩陣Table 1 Confusion matrix of RFI detection

(1)真陽(yáng)率(True positive rate, TPR):RFI被正確檢測(cè)的概率.

(2)假陽(yáng)率(False positive rate,FPR):非RFI數(shù)據(jù)中, 被標(biāo)記為RFI數(shù)據(jù)的概率.

(3)假陰率(False negative rate, FNR): RFI數(shù)據(jù)中, 被標(biāo)記為非RFI數(shù)據(jù)的概率.

(4)準(zhǔn)確率(Accuracy): 正確標(biāo)記的RFI數(shù)據(jù)和非RFI數(shù)據(jù)占全部數(shù)據(jù)的比例.

(5)精確率(Precision): 正確標(biāo)記為RFI的數(shù)據(jù)占全部標(biāo)記為RFI數(shù)據(jù)的比例.

SumThreshold算法在RFI檢測(cè)中具有較高的精度. 文獻(xiàn)[16]基于真陽(yáng)率和假陽(yáng)率對(duì)Sum-Threshold、SVD和VarThreshold等RFI檢測(cè)方法做了比較研究(詳見(jiàn)該文獻(xiàn)的圖8), 結(jié)果發(fā)現(xiàn)SumThreshold算法在RFI檢測(cè)中具有更高的精度.這主要是因?yàn)? SVD方法無(wú)法判斷數(shù)據(jù)幅度的平穩(wěn)增加是否是由RFI導(dǎo)致, 從而使它在檢測(cè)中易產(chǎn)生錯(cuò)誤.圖10為SumThreshold算法和VarThreshold算法的頻率-時(shí)間雙向RFI檢測(cè)結(jié)果對(duì)比, 由于RFI檢測(cè)策略不同, 使得SumThreshold算法和VarThreshold算法的RFI檢測(cè)結(jié)果不相同. 實(shí)際上, Var-Threshold算法的RFI檢測(cè)結(jié)果是SumThreshold算法RFI檢測(cè)結(jié)果的子集(圖10 (c)和(d)), SumThreshold算法更易檢測(cè)出各類(lèi)RFI數(shù)據(jù).

圖10 SumThreshold算法和VarThreshold算法的頻率-時(shí)間雙向RFI檢測(cè)結(jié)果對(duì)比. (a) SumThreshold算法(采用頻率-時(shí)間雙向檢測(cè)策略)的RFI檢測(cè)結(jié)果; (b) VarThreshold算法(采用頻率-時(shí)間雙向檢測(cè)策略)的RFI檢測(cè)結(jié)果; (c) SumThreshold算法檢測(cè)到但未被VarThreshold算法檢出的RFI; (d) VarThreshold算法檢測(cè)到但未被SumThreshold算法檢出的RFI. 該數(shù)據(jù)來(lái)源于FAST觀測(cè)[23].Fig.10 Comparison between the RFI detection results of SumThreshold and VarThreshold. (a) The RFI detection result using SumThreshold with frequency-time bidirection scheme; (b) the RFI detection result using VarThreshold with frequency-time bidirection scheme; (c) the RFI detected by SumThreshold method, but not by VarThreshold method; (d) the RFI detected by VarThreshold, but not by SumThreshold method. The experimental data are of a FAST observation[23].

基于2.3節(jié)RFI檢測(cè)的示例數(shù)據(jù), SumThreshold算法和VarThreshold算法的RFI檢測(cè)對(duì)比結(jié)果如圖11. 通過(guò)定量對(duì)比發(fā)現(xiàn), VarThreshold算法更容易出現(xiàn)漏標(biāo)記的現(xiàn)象. 同時(shí)文獻(xiàn)[12]指出Sum-Threshold算法假陽(yáng)率的定量化計(jì)算更為簡(jiǎn)單, 使得它在大規(guī)模觀測(cè)數(shù)據(jù)中檢測(cè)RFI時(shí)更為高效.

圖11 VarThreshold算法和SumThreshold算法RFI檢測(cè)對(duì)比圖. (a)一個(gè)觀測(cè)數(shù)據(jù)的截面; (b) VarThreshold算法RFI檢測(cè)結(jié)果; (c)SumThreshold算法RFI檢測(cè)結(jié)果.Fig.11 Comparison between the RFI detection results of VarThreshold and SumThreshold. (a) An intersecting surface of an observation; (b) the intersecting surface of the RFI detection result using VarThreshold method; (c) the intersecting surface of the RFI detection result using SumThreshold method.

3 SumThreshold算法優(yōu)化

3.1 形態(tài)學(xué)算子

在有些情況下, 干擾源的接收功率會(huì)隨著時(shí)間和頻率產(chǎn)生變化. 這使得即使連續(xù)地接收同一個(gè)干擾源的數(shù)據(jù), 也會(huì)因接收功率的變化, 導(dǎo)致閾值檢測(cè)方法無(wú)法在整個(gè)范圍內(nèi)檢測(cè)到干擾源. 圖12是SumThreshold算法未做膨脹操作的RFI檢測(cè)結(jié)果. 該檢測(cè)結(jié)果中存在大量干擾源接收功率變化造成的帶狀干擾斷裂. 對(duì)此, 可利用形態(tài)學(xué)方法將高聚集RFI數(shù)據(jù)周?chē)臄?shù)據(jù)進(jìn)一步做RFI標(biāo)記, 消除斷裂.

圖12 SumThreshold算法未做膨脹操作的RFI檢測(cè)結(jié)果. 該檢測(cè)結(jié)果中存在大量干擾源接收功率變化造成的帶狀干擾斷裂. 該數(shù)據(jù)來(lái)源于FAST觀測(cè)[23].Fig.12 The RFI detection results using the SumThresold without morphology-based flagging. The figure shows negative effects from the unstableness of received power from interfering sources. The experimental data are of a FAST observation[23].

文 獻(xiàn)[19]在 對(duì)LOFAR (Low Frequency Array)觀測(cè)數(shù)據(jù)進(jìn)行RFI檢測(cè)的研究中發(fā)現(xiàn), 普通的膨脹操作對(duì)尖銳的射頻特征較為敏感, 易將非射頻干擾數(shù)據(jù)誤判為RFI; 同時(shí), 通常的膨脹操作對(duì)平滑的射頻特征不敏感, 使得RFI數(shù)據(jù)存在漏標(biāo)現(xiàn)象. 因此, 傳統(tǒng)膨脹方法無(wú)法有效解決RFI檢測(cè)中的斷裂或邊緣漏檢測(cè)問(wèn)題. 因此, 文獻(xiàn)[19]引入新的膨脹算子: 尺度不變秩(Scale Invariant Rank,SIR)算子. 在使用SumThreshold算法進(jìn)行RFI檢測(cè)后, 進(jìn)一步基于SIR算子對(duì)檢測(cè)結(jié)果做后處理, 在時(shí)域、頻率或是兩者結(jié)合的方向上選取若干個(gè)大小不一的滑動(dòng)窗口逐行掃描, 計(jì)算滑動(dòng)窗口內(nèi)的置信度. 置信度是判斷滑動(dòng)窗口內(nèi)樣本是否應(yīng)該進(jìn)一步被標(biāo)記為RFI的指標(biāo). 當(dāng)置信度大于指定閾值時(shí), 將當(dāng)前滑動(dòng)窗口內(nèi)的樣本標(biāo)記為RFI. 反之, 如果置信度小于給定的閾值, 則當(dāng)前窗口內(nèi)的樣本不會(huì)被標(biāo)記為RFI. 置信度更新的準(zhǔn)則為: 從左到右對(duì)滑動(dòng)窗口內(nèi)樣本進(jìn)行遍歷, 當(dāng)樣本已被標(biāo)記為RFI時(shí), 置信度值加1, 反之, 置信度值不變.

SIR算子的引入進(jìn)一步提高了RFI檢測(cè)結(jié)果的精度. 基于高斯頻率分布的寬帶RFI特征模擬如圖13所示.圖13中當(dāng)圖像中加入噪聲,影響了RFI的可辨識(shí)度; 這時(shí)在噪聲數(shù)據(jù)上運(yùn)行SumThreshold算法會(huì)產(chǎn)生RFI標(biāo)識(shí)的漏標(biāo), RFI檢測(cè)準(zhǔn)確率為41.56%. 在SumThreshold算法檢測(cè)基礎(chǔ)上, 引入SIR算子, 通過(guò)圖13 (c)和(d)對(duì)比可以看出, SIR算子提高了檢測(cè)的精度, RFI檢測(cè)準(zhǔn)確率提高為59.53%. 但是不足在于, SIR運(yùn)算顯著地增加了算法的運(yùn)行時(shí)間.

圖13 基于高斯頻率分布的寬帶RFI特征模擬圖. (a)單獨(dú)RFI圖像; (b)加入噪聲的RFI圖像; (c)基于SumThreshold的RFI檢測(cè)結(jié)果; (d)基于SumThreshold + SIR操作的RFI檢測(cè)結(jié)果. 基于文獻(xiàn)[19]設(shè)計(jì)完成本實(shí)驗(yàn).Fig.13 Simulation of a typical broadband RFI feature with Gaussian frequency profile. (a) Isolated RFI feature; (b) when noise is added, a part of the feature becomes undetectable; (c) flagged with the SumThreshold method; (d) with SIR operator applied.This experiment is designed based on Ref.[19].

3.2 多特征檢測(cè)

SumThreshold算法主要是根據(jù)時(shí)頻觀測(cè)數(shù)據(jù)中的信號(hào)強(qiáng)度進(jìn)行RFI檢測(cè). 實(shí)際上, 觀測(cè)數(shù)據(jù)具有多種特征, 例如, 時(shí)間累計(jì)圖和相位等, 這些特征中均包含數(shù)據(jù)的有用信息. RFI數(shù)據(jù)與非RFI數(shù)據(jù)由于性質(zhì)差異, 在不同特征上會(huì)表現(xiàn)出不同的數(shù)據(jù)特點(diǎn). 例如, 文獻(xiàn)[16]研究發(fā)現(xiàn), 未污染數(shù)據(jù)中相位接近零旋轉(zhuǎn), 而RFI數(shù)據(jù)具有偏移為零的相位.因此, 相位包含有價(jià)值的RFI檢測(cè)信息, 針對(duì)多特征進(jìn)行基于SumThreshold的RFI檢測(cè)可提高算法精度.

3.3 基于圖形處理器(Graphics processing unit, GPU)的算法移植

隨著圖形處理元件的發(fā)展, 圖形處理器的應(yīng)用越來(lái)越廣泛. GPU計(jì)算相比于中央處理器(Central processing unit, CPU)具有并行度高、內(nèi)存帶寬高和運(yùn)行速度快等優(yōu)勢(shì). 因此, 已有學(xué)者[15]將SumThreshold算法移植到GPU平臺(tái). 由于體系結(jié)構(gòu)之間的差異, 在GPU平臺(tái)中對(duì)原始算法進(jìn)行了修改, 實(shí)現(xiàn)了基于GPU的小粒度并行處理. 基于GPU的改進(jìn)算法在速度和精度上均具有優(yōu)良的性能, 可實(shí)現(xiàn)RFI的在線檢測(cè).

4 總結(jié)

SumThreshold算法作為典型的閾值類(lèi)RFI檢測(cè)方法, 在射電數(shù)據(jù)處理中已有廣泛的應(yīng)用. 本文對(duì)SumThreshold算法進(jìn)行了詳細(xì)闡述, 介紹了算法的原理和流程, 分析了算法的關(guān)鍵步驟和設(shè)計(jì)實(shí)現(xiàn). 通過(guò)檢測(cè)示例, 直觀、詳盡地討論了SumThreshold算法進(jìn)行RFI檢測(cè)的過(guò)程. 性能分析表明, 此算法可有效檢測(cè)各類(lèi)型RFI, 滿足在線分析等高性能應(yīng)用需求. 針對(duì)算法的優(yōu)化, 本文從形態(tài)學(xué)算子、多特征選取、算法移植等角度分別做了總結(jié)和論述.

在未來(lái)的工作中, 將根據(jù)關(guān)鍵科學(xué)數(shù)據(jù)處理的需求進(jìn)行SumThreshold算法的改進(jìn)研究, 尤其是FAST數(shù)據(jù)的RFI檢測(cè). FAST作為世界上口徑最大和最靈敏的單口徑射電望遠(yuǎn)鏡, 它為眾多科學(xué)發(fā)現(xiàn)提供了前所未有的機(jī)遇. 在國(guó)家重大需求方面,FAST具有重要的應(yīng)用價(jià)值[25]. 因此, 有必要研究SumThreshold算法在FAST數(shù)據(jù)RFI 檢測(cè)方面的應(yīng)用.

致謝感謝審稿人對(duì)文章提出的寶貴建議, 使得文章的質(zhì)量有了顯著的提高.