林天琪 吳亞軍 朱人杰 徐志駿

(1中國科學(xué)院上海天文臺上海200030)

(2中國科學(xué)院大學(xué)北京100049)

(3中國科學(xué)院射電天文重點實驗室南京210023)

1 引言

近年來射電天文快速發(fā)展,一方面,國內(nèi)上海天馬65 m射電望遠鏡、500 m口徑球面射電望遠鏡(FAST)相繼投入運行,新疆奇臺110 m射電望遠鏡、云南景東120 m脈沖星射電望遠鏡等正在預(yù)研.另一方面,國際上平方公里陣(SKA)項目開始建造,其他大量的望遠鏡投入使用和升級,如:阿塔卡馬大型毫米波/亞毫米波陣列(Atacama Large Millimeter/Submillimeter Array,ALMA)、下一代的甚大天線陣(Next-Generation VLA,ngVLA)、國際低頻陣列(Low-Frequency Array,LOFAR)、默奇森寬場陣列(Murchinson Widefield Array,MWA)、加拿大CHIME(Canadian Hydrogen Intensity Mapping Experiment)天文望遠鏡等.與此同時超寬帶致冷接收系統(tǒng)的研究也取得進展,這些技術(shù)將射電觀測的靈敏度推向前所未有的高度,在此背景下射頻干擾(RFI)成為一個日益突出的問題.

對RFI緩解技術(shù)的方法可分為3類[1]:第1類是防止RFI信號進入天文數(shù)據(jù).在《無線電規(guī)則》和《中華人民共和國無線電頻率劃分規(guī)定》文件中,都通過對頻段和發(fā)射功率的相應(yīng)要求來保護射電天文業(yè)務(wù).在望遠鏡選址工作中,也會盡可能挑選遠離人類活動的區(qū)域,并建立無線電寧靜區(qū)域[2].第2類是實時去除數(shù)據(jù)中的RFI信號.實時去除RFI信號,可以在數(shù)據(jù)接收后對數(shù)據(jù)進行基于軟件和硬件的實時處理操作.WSRT(Westerbork Synthesis Radio Telescope)利用FPGA(Field Programmable Gate Array)執(zhí)行初始頻域濾波,然后利用累積求和對接收到的數(shù)據(jù)進行閾值處理[3].在相控陣中利用輔助天線、自適應(yīng)調(diào)零等方法也可以很好地消除RFI信號[4].第3類是在相關(guān)和數(shù)據(jù)集成或數(shù)據(jù)緩沖之后去除RFI.現(xiàn)今對射頻干擾后相關(guān)消除的技術(shù)繁多,文獻[5–6]提出了利用形態(tài)運算符、SumThreshold方法、背景擬合技術(shù)相結(jié)合,組合成一個通用方法來處理RFI信號.該方法成功地用于幾個干涉望遠鏡,包括LOFAR、WSRT、VLA(Very Large Array)、GMRT(Giant Metrewave Radio Telescope)、ATCA(Australia Telescope Compact Array)、Arecibo、MWA以及單碟望遠鏡Parkes.文獻[7–9]中提出了一種基于譜峰度的RFI提取算法,通過判斷數(shù)據(jù)的譜峰度(SK)來識別RFI信號并標記.該算法利用得到的SK估計值的期望統(tǒng)計方差,定義了RFI檢測的閾值.

上述幾種方法能夠很好地根據(jù)數(shù)據(jù)特性對數(shù)據(jù)是否有干擾進行標記,但不能實現(xiàn)數(shù)據(jù)干擾的消減.而在文獻[10]中利用了獨立成分分析法來消除射頻干擾信號,將射頻干擾信號從信號中分離.獨立成分分析能夠很好地消除干擾信號,但需要多通道的數(shù)據(jù)進行分析.本文提出了基于2維小波變換的方法標記和消除射頻干擾,對信號中存在的干擾進行標記,并能夠?qū)⑿盘栔械纳漕l干擾信號去除.將數(shù)據(jù)的時間頻率序列作為原始數(shù)據(jù),用小波變換中的系數(shù)作為干擾存在的判據(jù)對數(shù)據(jù)進行標記,利用絕對中位差作為閾值將干擾信號分辨出來,通過改變變換后的小波系數(shù)對數(shù)據(jù)進行處理.實驗證明,該方法提高了數(shù)據(jù)的信噪比,能夠很好地標記干擾信號并進行消除.

2 理論基礎(chǔ)

2.1 射頻干擾信號分析

對于接收到的天文信號,可以建立一個這樣的模型:

其中S(n)為所需的天文信號,W(n)是系統(tǒng)接受的噪聲信號,I(n)為RFI信號.對于接受到的信號X(n),n為采樣序號.在沒有經(jīng)過處理前,所需的天文信號會淹沒在噪聲信號中,對于RFI信號的消除可以提高信噪比,使得天文信號在后續(xù)的處理中受到RFI信號的影響更小.

射電望遠鏡接收信號時可能會接收到發(fā)生在局部頻率范圍內(nèi),但在時間上具有連續(xù)性的信號,這類信號就是窄帶RFI.移動通信、GPS衛(wèi)星導(dǎo)航、電視廣播等都屬于窄帶RFI,其數(shù)學(xué)表達式為:

其中INB為窄帶RFI信號,j是虛數(shù)單位,A、f、θ分別為干擾信號的幅值、頻率和相位,這類信號在經(jīng)過傅里葉變換之后能夠清楚地分辨出來.在數(shù)據(jù)接收的過程中也會接收到一些寬帶的干擾信號.其中存在由于雷擊、電子設(shè)備短路等情況造成的瞬時寬帶RFI,也有輸電電纜、以太網(wǎng)電纜等持續(xù)的寬帶RFI信號.其數(shù)學(xué)表示為:

其中,式中IBB為寬帶RFI信號,Al、fl、θl分別為第l個干擾分量的幅值、頻率和相位,L為RFI源的個數(shù).

在實際的觀測中,望遠鏡接收到的信號可以用時間頻率序列來表示,RFI信號在圖像中也能夠很好地辨別出來.有短時間發(fā)生突變的較小干擾信號,也有貫穿整個圖像的寬帶持續(xù)信號.圖1(a)中豎線為寬帶脈沖RFI信號,橫線為窄帶RFI信號,圖1(b)中也存在寬帶持續(xù)時間信號,RFI往往并不是單一形式存在的,天文觀測中存在著復(fù)雜的電磁環(huán)境,RFI信號交疊存在,使得去除信號更為復(fù)雜.

圖1 寬帶及窄帶RFI信號的時間頻率序列.圖(a)中橫線為窄帶RFI信號,豎線為寬帶脈沖RFI信號,圖(b)是寬帶持續(xù)RFI信號.Fig.1 Time-frequency sequence of wideband and narrowband RFI signals.The horizontal line in panel(a)is the narrowband RFI signal,the vertical line in panel(a)is the broadband pulsed RFI signal,and panel(b)is the broadband continuous RFI signal.

2.2 2維離散小波變換

在處理時間頻率序列時,可以將其看作是一個2維圖像進行處理.由于小波變換能夠覆蓋整個頻域并且具有變焦特性,通過選擇合適的濾波器可以極大減小不同特征間的相關(guān)性,在不同頻段中可以利用不同的窗口對數(shù)據(jù)進行處理.在文獻[11]中,也利用了一維小波變換進行RFI信號消除.小波變換在圖像處理中得到了廣泛的應(yīng)用,因此選用2維離散小波變換對時間頻率圖像進行處理.

小波變換是一種信號的時間-尺度分析方法,是近代信號分析和處理的一種重要手段.小波變換在時頻兩域都具有表征信號局部特征的能力,是一種窗口大小固定不變但形狀、時間窗和頻率窗都可以改變的時頻局部分析方法.通過小波變換,對信號進行多尺度的細化分析,能夠有效地從信號中提取信息.

由于數(shù)據(jù)為2維信號f(x,y),因此采用2維離散小波變換對信號進行分解.2維離散小波變換需要一個2維尺度函數(shù)和3個2維小波ψH(x,y)、ψV(x,y)、ψD(x,y).ψH(x,y)、ψV(x,y)、ψD(x,y)分別表示信號沿行(H)、列(V)和對角線(D)方向的變化[12–13].

對于給定的尺度函數(shù)φr,p,q(x,y)和平移基函數(shù)(x,y):

其中,r為尺度參數(shù),p、q分別為沿著x、y方向上的平移參數(shù).對于像素為P×Q大小的圖像,經(jīng)過小波變換分解后的低頻系數(shù)為:

其中,wφ(r0,p,q)為分解后的低頻系數(shù),r0是分解后第0層的系數(shù).經(jīng)過小波變換分解后的高頻系數(shù)為:

2.3 RFI信號在小波系數(shù)中的體現(xiàn)

將數(shù)據(jù)進行小波變換后,可以得到低頻和高頻的小波系數(shù),由于RFI存在的形式多種多樣,在系數(shù)中體現(xiàn)也不盡相同.

對于持續(xù)寬帶RFI信號,經(jīng)過小波變換后RFI能量會更多地體現(xiàn)在小波變換的低頻部分.在圖像中加入高斯信號作為原始信號,再加入一定寬度的信號作為寬帶RFI信號.經(jīng)過小波變換,得到分解后的低頻和高頻系數(shù)如圖2所示.由于寬帶信號分布映射到時間頻率序列圖像中占用的面積較大,所以信號分解的小波系數(shù)能量會集中在信號的低頻部分,但是在加入的信號邊緣處,也會在小波分解后的高頻系數(shù)中存在.

圖2 窄帶RFI在小波變換中的體現(xiàn).上圖為原始信號,中圖為小波分解后的低頻信號,下圖為小波分解后的第3層高頻信號.Fig.2 The embodiment of narrowband RFI in wavelet transform.The top panel is the original signal,the middle panel is the low frequency signal after wavelet decomposition,and the bottom panel is the third layer of high frequency signal after wavelet decomposition.

而對于RFI信號中的寬帶脈沖信號、窄帶信號以及只有一個突變值的點頻瞬時信號,在時間頻率序列中往往表現(xiàn)為一條直線和一個很亮的點.它們的信號不一定會出現(xiàn)在低頻信號中,但會因為具有突變性而在高頻系數(shù)中得以體現(xiàn).

在高斯圖像中加入窄帶信號和帶有突變值的點信號,經(jīng)過小波變換后得到低頻和高頻系數(shù),如圖3所示.當窄帶信號較強時,信號也可以在低頻系數(shù)中有所體現(xiàn),但對于信號較弱或者不連續(xù)時,會在高頻系數(shù)中體現(xiàn)得更加明顯.

本文采用絕對中值濾波的方式對于小波變換后的系數(shù)進行處理.由于天文信號大多淹沒在高斯噪聲中,所以可以近似看成高斯信號,絕對中位差(Median absolute deviation,MAD)作為檢測單變量樣本中樣本差異性的穩(wěn)健度量,相較于利用均值和標準差,MAD能夠更好地判別出信號中的異常值,可以大大減少其對于數(shù)據(jù)集的影響.在文獻[14–15]中,也利用了MAD方法對RFI信號進行去除.對于給定的信號X(x1,x2,···,xK),其中K為樣本個數(shù),其MAD值如下:

當分布為正態(tài)分布時,使用σr=1.4826×MAD作標準差σ的一致估計量.相較于樣本差和方差更不容易受到異常值的影響.

3 數(shù)據(jù)處理分析

3.1 RFI信號標記

由于干擾信號在低頻中大多能夠體現(xiàn)出來,對于RFI信號的標記,我們利用小波變換得到的低頻系數(shù)做閾值處理.為了確保標記的精度,應(yīng)盡量選擇較小的小波閾值層數(shù),將小波基函數(shù)與信號值進行卷積,通過處理標記低頻系數(shù)的閾值來標記RFI信號.

圖3 寬帶RFI在小波變換中的體現(xiàn).上圖為原始信號,中圖為小波分解后的低頻信號,下圖為小波分解后的第3層高頻信號.Fig.3 The embodiment of broadband RFI in wavelet transform.The top panel is the original signal,the middle panel is the low frequency signal after wavelet decomposition,and the bottom panel is the third layer of high frequency signal after wavelet decomposition.

利用時間頻率序列作為信號原始圖像,對信號進行處理得到小波系數(shù).若原始數(shù)據(jù)波動較大,可以先對數(shù)據(jù)做平滑濾波.處理后的數(shù)據(jù)如圖4所示,圖中的橫線為計算得到的閾值.從圖中能夠明顯看到有干擾的部分模值相對較大,其中超過閾值的部分將被標注出來.數(shù)據(jù)標注后,將小波系數(shù)還原,選擇處理前后變化值較大的數(shù)據(jù)進行標記.

快速射電暴(FRB)是一種高能量、持續(xù)時間為毫秒的無線電脈沖.不考慮傳播中的閃爍和散射,FRB是一個持續(xù)時間短的寬帶結(jié)構(gòu),與脈沖星中觀察到的單脈沖類似,當FRB穿過電離層時,信號會產(chǎn)生與頻率相關(guān)的時延,傳播信號的長波比短波要慢上許多,經(jīng)過消色散的處理之后,可以通過不同頻率上數(shù)據(jù)的疊加來識別FRB的存在.RFI也會干擾FRB的識別.圖5是一個帶有RFI的FRB數(shù)據(jù),數(shù)據(jù)有336個通道,觀測波段為1120–1465 MHz,采樣間隔為0.0012 s.對數(shù)據(jù)進行標記后,通過填充隨機數(shù)的方式,對標記的干擾信號進行處理,得到進行數(shù)據(jù)標記處理和消除色散后的結(jié)果.圖中可以看出信噪比由標記之前的5.06變?yōu)?.60.

圖4 時間頻率序列及其系數(shù)閾值處理.上圖為原始信號,下圖為小波低頻系數(shù)與閾值的比較圖,圖中橫線為閾值.Fig.4 Time-frequency sequence and its coefficient threshold processing.The top panel is the original signal,the bottom panel is the comparison diagram of wavelet low-frequency coefficients and the threshold,and the horizontal line in the figure is the threshold.

3.2 RFI信號去除

在將信號圖像進行2維小波變換后,干擾信號具有更大的模值.在進行RFI信號去除的工作中,期望能在原始數(shù)據(jù)盡可能不變的情況下更好地消除RFI干擾信號.為了盡量保留更多的原始數(shù)據(jù),對于低頻系數(shù),先將得到的超出閾值的系數(shù)進行低通濾波,再與原始系數(shù)相減,使得更多的信號被留存下來.

干擾信號的去除過程如圖6所示,具體的操作步驟如下:

(1)對要處理的信號圖像進行多層小波變換后,得到各層的小波系數(shù);

(2)對于低頻系數(shù),需要對得到的小波系數(shù)進行低通濾波,去除有干擾的信號,留下變化較大的信號.圖7為系數(shù)處理的過程,其中標記為藍色的線為小波分解后的系數(shù),紅色的線為濾波后的信號,黃色的線是小波系數(shù)與濾波后結(jié)果的差值;

(3)將處理后的小波系數(shù)進行閾值處理,選擇MAD閾值λ,將高于閾值的系數(shù)置零,得到估計的小波系數(shù);

(4)處理后的小波系數(shù)重構(gòu),得到去除干擾信號的圖像.

圖6 RFI信號去除過程Fig.6 RFI signal removal process

圖7 小波系數(shù)處理.藍色為原始信號,紅色為濾波后的信號,黃色為最終結(jié)果.Fig.7 Wavelet coefficient processing.Blue is the original signal,red is the filtered signal,and yellow is the final result.

在去除時間頻率序列的RFI信號時,首先對數(shù)據(jù)進行小波變換,將標記的低頻系數(shù)進行平滑處理.然后對標記的高頻系數(shù)進行閾值處理,最后將信號進行還原.還原信號的結(jié)果如圖8.圖8是一個同時帶有FRB信號和噪聲信號的數(shù)據(jù),數(shù)據(jù)有221個通道,觀測波段為2189.5–2300 MHz,采樣間隔為0.002 s.可以看到經(jīng)過處理后噪聲數(shù)據(jù)被消減,并且保留了FRB信號,信噪比從7.56變?yōu)?.11.

圖8 RFI信號去除結(jié)果.上圖是處理前的圖像,下圖是處理后的圖像,MAX為數(shù)據(jù)的最高信噪比.Fig.8 The result of RFI signal removal.The top panel is the image before processing,the bottom panel is the image after processing,and MAX is the maximum signal-to-noise ratio of the data.

4 總結(jié)與展望

本文提出了用小波變換來標記和去除時間頻率序列中的干擾信號,數(shù)據(jù)通過小波變換之后,利用低頻系數(shù)標記出干擾信號,并對小波系數(shù)中的高低頻系數(shù)分別處理,利用絕對中位差作為閾值來消除干擾信號.通過對仿真數(shù)據(jù)和實際觀測的數(shù)據(jù)進行測試,證實該方法能夠提高數(shù)據(jù)的信噪比.由于在信號判別過程中使用的是小波變換的低頻系數(shù),如果分解的層數(shù)過多,很有可能會擴大數(shù)據(jù)的標記范圍.而對于一些強度不高的信號點,有可能未被標記出來.在干擾數(shù)據(jù)與信號數(shù)據(jù)重疊的時候,進行隨機數(shù)據(jù)替換可能降低信號的信噪比.射頻干擾信號紛繁復(fù)雜,利用閾值處理的結(jié)果仍然會有誤判和少判的情況存在,后續(xù)可以繼續(xù)優(yōu)化閾值處理方法,如利用系數(shù)間相關(guān)性去除干擾信號.