基于小波變換的射頻干擾消除方法的研究*

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

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

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

(3中國科學(xué)院射電天文重點(diǎn)實(shí)驗(yàn)室南京210023)

1 引言

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

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

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

2 理論基礎(chǔ)

2.1 射頻干擾信號(hào)分析

對于接收到的天文信號(hào),可以建立一個(gè)這樣的模型:

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

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

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

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

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

圖1 寬帶及窄帶RFI信號(hào)的時(shí)間頻率序列.圖(a)中橫線為窄帶RFI信號(hào),豎線為寬帶脈沖RFI信號(hào),圖(b)是寬帶持續(xù)RFI信號(hào).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維離散小波變換

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

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

由于數(shù)據(jù)為2維信號(hào)f(x,y),因此采用2維離散小波變換對信號(hào)進(jìn)行分解.2維離散小波變換需要一個(gè)2維尺度函數(shù)和3個(gè)2維小波ψH(x,y)、ψV(x,y)、ψD(x,y).ψH(x,y)、ψV(x,y)、ψD(x,y)分別表示信號(hào)沿行(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信號(hào)在小波系數(shù)中的體現(xiàn)

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

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

圖2 窄帶RFI在小波變換中的體現(xiàn).上圖為原始信號(hào),中圖為小波分解后的低頻信號(hào),下圖為小波分解后的第3層高頻信號(hào).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信號(hào)中的寬帶脈沖信號(hào)、窄帶信號(hào)以及只有一個(gè)突變值的點(diǎn)頻瞬時(shí)信號(hào),在時(shí)間頻率序列中往往表現(xiàn)為一條直線和一個(gè)很亮的點(diǎn).它們的信號(hào)不一定會(huì)出現(xiàn)在低頻信號(hào)中,但會(huì)因?yàn)榫哂型蛔冃远诟哳l系數(shù)中得以體現(xiàn).

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

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

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

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

3.1 RFI信號(hào)標(biāo)記

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

圖3 寬帶RFI在小波變換中的體現(xiàn).上圖為原始信號(hào),中圖為小波分解后的低頻信號(hào),下圖為小波分解后的第3層高頻信號(hào).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í)間頻率序列作為信號(hào)原始圖像,對信號(hào)進(jìn)行處理得到小波系數(shù).若原始數(shù)據(jù)波動(dòng)較大,可以先對數(shù)據(jù)做平滑濾波.處理后的數(shù)據(jù)如圖4所示,圖中的橫線為計(jì)算得到的閾值.從圖中能夠明顯看到有干擾的部分模值相對較大,其中超過閾值的部分將被標(biāo)注出來.數(shù)據(jù)標(biāo)注后,將小波系數(shù)還原,選擇處理前后變化值較大的數(shù)據(jù)進(jìn)行標(biāo)記.

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

圖4 時(shí)間頻率序列及其系數(shù)閾值處理.上圖為原始信號(hào),下圖為小波低頻系數(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信號(hào)去除

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

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

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

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

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

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

圖6 RFI信號(hào)去除過程Fig.6 RFI signal removal process

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

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

圖8 RFI信號(hào)去除結(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é)與展望

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