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

(1 華南師范大學計算機學院 廣州 510631)

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

1 引言

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

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

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

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

2 SumThreshold算法描述

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

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

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

2.1 基線矯正

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

圖2 基線對閾值類方法的影響. 這是FAST的一個時間-頻率觀測圖像, 觀測時間為0.05 s, 頻率范圍是1000–1500 MHz[23]. (a)原始觀測數(shù)據(jù);(b)去除基線效應后的結果. 原始觀測數(shù)據(jù)中不同頻帶像素讀數(shù)浮動范圍變化非常大, 對RFI檢測中的閾值設定造成困擾.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.

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

2.2 RFI檢測方向指定

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

圖3 FAST觀測數(shù)據(jù). 觀測時間為0.05 s, 頻率范圍是1000–1500 MHz. 由于存在取值較大的RFI, 導致圖像細節(jié)不明顯.該數(shù)據(jù)來源于FAST觀測[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].

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

圖4 SumThreshold (ST)算法的不同實施策略效果比較. (a) ST算法沿時間維度的RFI檢測結果; (b) ST算法沿頻率維度的RFI檢測結果;(c)在圖(a)中被檢測到但未在圖(b)中檢測到的RFI; (d)在圖(b)中被檢測到但未在圖(a)中檢測到的RFI. 圖中白色點表示RFI數(shù)據(jù). 由于強干擾被兩種策略均正確檢測, 所以在圖(c)和(d)中均未標記, 且它們比圖(a)和圖(b)整體上顯得暗一些. 該數(shù)據(jù)來源于FAST觀測[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檢測結果. (a)先沿著時間方向再沿著頻率方向的RFI檢測結果(簡稱為時間-頻率雙向檢測); (b)先沿著頻率方向再沿著時間方向的RFI檢測結果(簡稱為頻率-時間雙向檢測); (c)將圖4 (a)和(b)單一方向檢測結果合并后的RFI標記結果. 該數(shù)據(jù)來源于FAST觀測[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 基于不同雙向檢測策略的SumThreshold算法效果對比. (a)被時間-頻率雙向方式檢測到但未被頻率-時間雙向檢測方式發(fā)現(xiàn)的RFI; (b)被頻率-時間雙向檢測方式發(fā)現(xiàn)但未被時間-頻率雙向檢測方式發(fā)現(xiàn)的RFI. 由于強干擾被兩種雙向策略均正確檢測, 所以它們在這個實驗中均未標記, 且該圖比圖5整體上暗一些. 該數(shù)據(jù)來源于FAST觀測[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算法對RFI進行檢測

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

2.3.1 RFI數(shù)據(jù)標記準則

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

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

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

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

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

2.3.2 SumThreshold算法設計

對于某個給定寬度為Mn的檢測窗口, Sum-Threshold檢測的過程請見算法1:

算法1:

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

輸出: 標記集合{f0,f1,··· ,fL-1}

(1)獲取當前窗口大小Mn、閾值χn、初始化變量i←0、未檢測為RFI的數(shù)據(jù)之和Sum←0、Count←0, 本次迭代生成的標記集合{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), 計 算 從 第1個 數(shù) 據(jù) 開 始, 長為Mn的窗口內(nèi)尚未被檢測為RFI的像素的均值.步驟(10)–(14)判斷當前窗口內(nèi)像素的均值是否大于閾值, 如果判斷結果為“是”, 則將當前窗口內(nèi)所有尚未被檢測為射頻干擾的像素標記為RFI. 步驟(15)–(18)的意思是, 如果檢測窗口外右側像素在之前的迭代中未被判斷為RFI數(shù)據(jù), 則將此數(shù)據(jù)移入窗口, 窗口內(nèi)數(shù)據(jù)之和以及數(shù)據(jù)個數(shù)增加. 步驟(19)–(22)的意思是如果當前窗口內(nèi)左側像素在之前的迭代中未被判斷為RFI數(shù)據(jù), 則將其移出窗口, 窗口內(nèi)數(shù)據(jù)之和以及數(shù)據(jù)個數(shù)相應減少. 循環(huán)執(zhí)行步驟(10)–(22), 直至所有數(shù)據(jù)被處理, 更新標記集合, 本次迭代結束.

2.3.3 基于SumThreshold算法的RFI檢測示例

圖8顯示了一個觀測數(shù)據(jù)的截面及其RFI標記結果的示例. 該數(shù)據(jù)是從某個頻帶中截取得到. 因此, 圖8 (a)的橫坐標和縱坐標分別表示時間和流量. 經(jīng)過基線擬合和剔除后, 各像素的值分別為:{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算法進行RFI檢測后, 值為{12,14, 16, 15, 17, 20, 24, 26, 22, 18, 14, 13, 11}的像素被標記為RFI (圖8 (b)中斜線標記的部分).

圖8 一個觀測數(shù)據(jù)及其RFI檢測結果的截面. (a)一個觀測數(shù)據(jù)的截面; (b)左圖數(shù)據(jù)的RFI檢測結果.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檢測過程如圖9所示. 在該實驗中, SumThreshold算法經(jīng)過了3次迭代, 第1次迭代滑動窗口寬度是M1= 1, 閾值為χ1= 20. 迭代輸入數(shù)據(jù)為圖8 (a)中的數(shù)據(jù), 如圖9 (a)所示. 因為窗口大小為1, 隨著窗口滑動, 每一個數(shù)據(jù)均與閾值進行比較, 值為24、26和22的數(shù)據(jù)大于閾值. 因此, 相應像素被標記為RFI (圖9 (b)). 第2次迭代的滑動窗口寬度為M2= 2, 閾值是χ2= 13.3333.為了防止已標記為RFI的數(shù)據(jù)值過大導致周圍的像素被錯誤檢測為射頻干擾, 在第2次迭代檢測之前將第1輪迭代時被標記為RFI的像素數(shù)值改為當前閾值13.3333 (圖9 (c)). 第2次迭代后新被標記為RFI的數(shù)據(jù)分別為14、16、15、17、20、18、14和13 (圖9 (d)). 對于值為13的數(shù)據(jù), 當其進入滑動窗口時, 窗口內(nèi)的數(shù)據(jù)值為14和13, 該窗口內(nèi)像素的均值大于閾值13.3333, 因此標記數(shù)據(jù)13為RFI. 如果將值為13的像素單獨與閾值相比,則由于小于閾值而不會被檢測到, 但是因為窗口內(nèi)像素的整體平均效應, 該像素被SumThreshold算法標記為RFI, 解決了傳統(tǒng)閾值法RFI由單個像素比較導致的漏檢問題. 第3次迭代中窗口寬度是M3= 4, 檢測閾值為χ3= 10.5180. 之前迭代中被標記為RFI的像素的值在本次迭代中置換為10.5180 (圖9 (e)). 通過與第2次迭代相同操作流程, 滑動窗口對數(shù)據(jù)處理后, 新被標記為RFI的數(shù)據(jù)是值為12和11的像素(圖9 (f)).

圖9 RFI檢測過程示例圖. 圖(a)、(c)、(e)分別是第1、2、3次迭代前的輸入數(shù)據(jù); 圖(b)、(d)、(f)分別是第1, 2, 3次迭代處理后的數(shù)據(jù). 詳述見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ù)據(jù)長度為L, 則理論上檢測窗口大小的集合為{1,2,··· ,L-1}. 每次迭代要求對每個窗口子序列大小的數(shù)據(jù)進行傳遞. SumThreshold算法的時間復雜度為O(L2). 這對于算法效率要求較高的應用很難滿足要求. 因此, 基于指數(shù)增長設定窗口大小, 則檢測窗口集合成為原始窗口集合的子集. 例如, 若將窗口大小設為[1,2,4,8,16,···],則算法時間復雜度為O(Llog2L). 與理論上的檢測窗口集合相比, 減少了尺寸為[3,5,6,7,···]的窗口.這些被減少的窗口所檢測的數(shù)據(jù)特征, 仍可能被約簡窗口集合檢測到. 因此, 上述對候選窗口減少的做法對算法精度不產(chǎn)生顯著影響.

在減少候選檢測窗口集合的基礎上, 設定數(shù)據(jù)集子集大小, 可進一步將算法復雜度降低. 在文獻[24]中, 通過限定數(shù)據(jù)子集小于1024, 算法效率明顯提高. 但是算法的精度會受影響, 這是因為非常微弱或者范圍大的特征將被忽略.

2.4.2 算法適用性

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

2.4.3 算法精度

SumThreshold算法對射電數(shù)據(jù)進行RFI標記,結果可分為4類: 真陽性(True positive, TP): 標記為RFI的數(shù)據(jù), 實際上也是RFI數(shù)據(jù); 假陽性(False positive, FP): 標記為RFI的數(shù)據(jù), 實際上不是RFI數(shù)據(jù); 真陰性(True negative, TN): 標記為非RFI的數(shù)據(jù), 實際上也是非RFI數(shù)據(jù); 假陰性(False negative, FN): 標記為非RFI的數(shù)據(jù), 實際上是RFI數(shù)據(jù). 這4類結果可構成混淆矩陣, 見表1.

基于表1, 射頻干擾檢測結果的常用精度度量指標有:

表1 RFI檢測混淆矩陣Table 1 Confusion matrix of RFI detection

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

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

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

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

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

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

圖10 SumThreshold算法和VarThreshold算法的頻率-時間雙向RFI檢測結果對比. (a) SumThreshold算法(采用頻率-時間雙向檢測策略)的RFI檢測結果; (b) VarThreshold算法(采用頻率-時間雙向檢測策略)的RFI檢測結果; (c) SumThreshold算法檢測到但未被VarThreshold算法檢出的RFI; (d) VarThreshold算法檢測到但未被SumThreshold算法檢出的RFI. 該數(shù)據(jù)來源于FAST觀測[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檢測的示例數(shù)據(jù), SumThreshold算法和VarThreshold算法的RFI檢測對比結果如圖11. 通過定量對比發(fā)現(xiàn), VarThreshold算法更容易出現(xiàn)漏標記的現(xiàn)象. 同時文獻[12]指出Sum-Threshold算法假陽率的定量化計算更為簡單, 使得它在大規(guī)模觀測數(shù)據(jù)中檢測RFI時更為高效.

圖11 VarThreshold算法和SumThreshold算法RFI檢測對比圖. (a)一個觀測數(shù)據(jù)的截面; (b) VarThreshold算法RFI檢測結果; (c)SumThreshold算法RFI檢測結果.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)學算子

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

圖12 SumThreshold算法未做膨脹操作的RFI檢測結果. 該檢測結果中存在大量干擾源接收功率變化造成的帶狀干擾斷裂. 該數(shù)據(jù)來源于FAST觀測[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].

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

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

圖13 基于高斯頻率分布的寬帶RFI特征模擬圖. (a)單獨RFI圖像; (b)加入噪聲的RFI圖像; (c)基于SumThreshold的RFI檢測結果; (d)基于SumThreshold + SIR操作的RFI檢測結果. 基于文獻[19]設計完成本實驗.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 多特征檢測

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

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

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

4 總結

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

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

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