馬 琴 裴 鑫 陳卯蒸 托乎提努爾 王兆軍 李 健段雪峰

(1新疆大學物理科學與技術學院 烏魯木齊 830046)

(2中國科學院新疆天文臺 烏魯木齊 830011)

(3中國科學院大學 北京 100049)

(4新疆微波技術重點實驗室 烏魯木齊 830011)

1 引言

隨著科技的發(fā)展,射電望遠鏡的靈敏度、觀測視場以及空間分辨率等諸多性能都有著比較明顯的提升.但與此同時,無線電業(yè)務領域的應用范圍也在不斷擴大,使得射電天文望遠鏡不可避免地受到非天文信號的影響,即射頻干擾(Radio Frequency Interference,RFI),其主要來源包括無線通信、無線網(wǎng)絡、雷達、廣播電視、衛(wèi)星和人們使用的各類電子設備[1].RFI不但使得觀測天文信號質量降低,而且將耗費大量的信號處理資源和存儲空間.因此如何避免或減少RFI被觀測系統(tǒng)接收和記錄,已經(jīng)成為射電望遠鏡在觀測中需要解決的突出問題.

傳統(tǒng)的單通道射電天文干擾消除方法[2–5]在信噪比較低時只能去除部分干擾,且容易造成觀測信號失真.基于自適應濾波器的干擾消除方法可以在信號和噪聲的統(tǒng)計特性未知的前提下,從噪聲中提取信號,當信號與噪聲的統(tǒng)計特性發(fā)生變化時,能調(diào)節(jié)自身的參數(shù)以適應不同狀況,從而實現(xiàn)最優(yōu)濾波的特性[6].它通過將射電天文信號自身特點與雷達及無線電通信等領域微弱信號處理方法相結合,提出了自適應干擾對消的策略.該方法借助一架低增益參考天線來實時監(jiān)測射頻干擾,然后由自適應濾波器跟蹤兩個輸入信號并去除RFI,保留所需的天文信號.

自適應干擾對消早在20世紀60年代就被引入,作為減少低頻系統(tǒng)中干擾噪聲的一種方法,成功應用于語音信號噪聲消除、心電圖周期噪聲信號抵消等領域[7].其最主要的部件是自適應濾波器,而自適應濾波器在普通濾波器的基礎上加入了調(diào)整濾波器系數(shù)的自適應算法,根據(jù)自適應算法優(yōu)化準則可選取不同的算法,但這些算法的計算量差異較大,故而在選擇算法時既要考慮性能也要考慮計算量的大小.自適應濾波器兩種常用的算法是最小均方(Least Mean Square,LMS)和遞推最小二乘法(Recursive Least Square,RLS)[8],RLS算法具有較快的收斂速度,但其運算量較大且具有發(fā)散性,LMS算法收斂速度一般,但收斂性好.本文將對這兩種方法進行分析和比較,重點采用LMS進行建模,利用仿真軟件對自適應濾波器進行設計和分析,以系統(tǒng)誤差和收斂性為評判標準,通過改變步長與階數(shù)對濾波效果進行優(yōu)化.

2 研究內(nèi)容

2.1 自適應干擾對消模型

自適應干擾對消系統(tǒng)是基于自適應濾波器原理的一種擴展[9],其原理如圖1所示,它由主信道(射電望遠鏡)和一個單獨的參考通道組成.主信道接收被RFI(無線電頻率干擾)污染的天文信號,即進入射電望遠鏡的天文信號s(n)和干擾iP(n),其中iP(n)也包含了白噪聲,n為時間序列.一架低增益參考天線被用來接收干擾iR(n),其中由于參考天線增益有限,因此在濾波器自適應時間尺度上的參考輸入中基本上沒有天文信號[10].

圖1 自適應干擾對消原理圖Fig.1 Schematic diagram of adaptive interference cancellation

參考信道iR(n)中的干擾與主信道iP(n)中的干擾相關,自適應濾波器工作的主要目的是將這種相關性作為時間的函數(shù).自適應算法通過比較前一時刻和當前的信息,并發(fā)送更新后的加權系數(shù)到數(shù)字濾波器,數(shù)字濾波器利用這些系數(shù)改變iR(n),產(chǎn)生與主信道干擾相似的輸出y(n);從一次輸入中減去y(n),得到系統(tǒng)輸出ε(n):

對于s(n)、iP(n)、iR(n)或它們之間的相互關系,不需要事先知道.自適應算法通過比較ε(n)和ε(n?1)找到新的系數(shù),使用LMS算法最小化總功率:

由于s(n)與主路和參考路的干擾不相關,交叉項為零,因此系統(tǒng)輸出的期望值(時間平均)為:

當濾波器對系數(shù)進行調(diào)整使E{ε2(n)}最小時,天文信號E{s2(n)}中的功率被吸收,因此E{[iP(n)?y(n)]2}達到功率最小,即:

由于天文信號是恒定的,使總輸出功率最小就使輸出的干擾功率最小,從而使輸出信噪比最大.此時的ε(n)就是s(n)的最佳估計,y(n)實際上是iP(n)在最小均方意義下的估計.

2.2 自適應濾波器結構

自適應濾波器分為有限沖激響應(Finite Impulse Response,FIR)和無限沖激響應(Infinite Impulse Response,IIR)這兩種結構[11].FIR濾波器由有限個脈沖響應函數(shù)離散值的濾波器組成,它在實際應用中最常用,可以實現(xiàn)較為嚴格的線性相移特性,并且可以保證濾波后波形不失真,但若要求頻域過渡帶快速衰減,就需要更多的階數(shù).IIR濾波器可以以較小的計算量獲得陡降的過渡帶,但其收斂慢、穩(wěn)定性差,且較難保證線性相移,在FIR濾波器中則不存在類似問題.FIR濾波器結構又可以分為3種結構類型:格型結構、橫向結構和對稱橫向結構,其中橫向結構因其形式簡單及易于實現(xiàn)最為常用,故本文采用FIR橫向濾波器結構進行設計,濾波器的結構如圖2所示.

濾波器的輸出為:

其中XT(n)為濾波器輸入矢量,即X(n)=[x(n)x(n?1)···x(n?N+1)]T,W(n)為權系數(shù)矢量(濾波器的沖激響應),即W(n)=[w0(n)w1(n)w2(n)w3(n)···w N?1(n)]T,N為濾波器的階數(shù)[12].

圖2 FIR濾波器結構Fig.2 The structure of FIR filter

2.3 自適應濾波器算法

本文采用LMS算法,因其具有計算量小、性能穩(wěn)健、易于實現(xiàn)、不依賴于模型等優(yōu)點,在實踐中被廣泛使用[13].

LMS算法是一種基于最小均方誤差準則的隨機梯度下降算法,核心思想就是利用均方誤差代替平方誤差,其通過調(diào)節(jié)權系數(shù)使得濾波器的均方誤差最小,誤差曲面的梯度為:

最陡下降法迭代計算權矢量的公式為:

式中μ是控制步長的參數(shù)稱為自適應增益常數(shù),是表征迭代速度快慢的物理量.W(n+1)和W(n)分別為迭代后和迭代前的系數(shù)值,從上式可以看出,自適應迭代下一時刻權系數(shù)矢量由上一時刻的權系數(shù)矢量加以誤差函數(shù)為收斂參數(shù)的輸入矢量得到.自適應濾波器收斂的條件是:

其中λmax是輸入信號自相關矩陣的最大特征值;當μ值越大自適應過程越快且時間越短,引起的失調(diào)越大,μ值越小,自適應過程越慢,系統(tǒng)越穩(wěn)定同時失調(diào)也越小,μ的選取必須在失調(diào)和收斂速度之間取得較好的折中,既要具有較快的收斂速度,又要使穩(wěn)態(tài)誤差最小.

算法步驟:

(1)初始化濾波器系數(shù),W(0)=0,或者可以根據(jù)先驗知識來確定初始權值.

(2)對時間序列每一時刻進行以下計算:

濾波:y(n)=XT(n)W(n),誤差:ε(n)=s(n)+iP(n)?iR(n),權向量更新:W(n+1)=W(n)+2με(n)X(n).

綜上所述,該方法的主要步驟為:(1)讀入主信道和參考信道的信號;(2)設置濾波器的階數(shù)和步長;(3)初始化自適應濾波算法的參數(shù);(4)自適應濾波處理;(5)濾波器權系數(shù)的更新.

3 仿真與分析

基于上述理論分析,由(5)式和(7)式可知,權系數(shù)矢量和濾波器輸入矢量乘積的和即為濾波器的估計輸出,加權系數(shù)矢量是通過步長、ε(n)以及與上一次迭代權值之和計算得到.根據(jù)圖1本文設計了不同階數(shù)和不同步長的自適應濾波器,對加干擾的正弦波進行濾波,輸入信號如圖3所示.

圖中,主通道中輸入目標信號、干擾以及白噪聲,參考信道輸入的只有干擾和白噪聲,通過LMS算法可以自適應調(diào)節(jié)線性組合器權系數(shù),主信道與參考信道內(nèi)的噪聲信號相互對消,輸出的信號即為所需的目標信號.

圖3 輸入信號Fig.3 The input signal

為獲取濾波器最優(yōu)參數(shù),本文首先分析了當階數(shù)一定時改變步長對濾波效果的影響,選取了階數(shù)為2、步長設置分別為0.0005、0.0002、0.00009的濾波器進行仿真;然后分析了當步長一定時改變階數(shù)對濾波效果的影響,選取了步長為0.0001、階數(shù)設置分別為2、6、12的濾波器進行仿真;同時用誤差對濾波性能進行了表征.如圖4和圖5分別是2階時不同步長下濾波后的信號和誤差,圖6和圖7分別是步長為0.0001時不同階數(shù)的濾波信號和誤差.

當階數(shù)一定,改變步長時,由圖4可知,計算結果隨步長增加收斂速度變快,步長越短的計算時間越長.濾波性能隨參數(shù)的改變在圖5中可以看到,誤差的收斂速度與濾波信號的收斂速度同步,但濾波效果隨步長增大而變差.

圖5 2階時不同步長的誤差(上:0.0005、中:0.0002、下:0.00009)Fig.5 The error of different step sizes when the order is 2(top panel:0.0005,middle panel:0.0002,bottom panel:0.00009)

圖5 續(xù)Fig.5 Continued

圖6 步長為0.0001時不同階數(shù)的濾波信號(上:12、中:6、下:2)Fig.6 Filtered signals of different orders when the step size is 0.0001(top panel:12,middle panel:6,bottom panel:2)

圖7 步長為0.0001時不同階數(shù)的誤差(上:12、中:6、下:2)Fig.7 The error of different orders when the step size is 0.0001(top panel:12,middle panel:6,bottom panel:2)

當步長一定改變階數(shù)時,如圖6所示,濾波信號收斂時間隨著階數(shù)的減小而增大,濾波信號的收斂速度隨著階數(shù)的增大而變快.濾波效果如圖7所示,隨著階數(shù)的減小,誤差也更接近0.

在進行多次仿真實驗后發(fā)現(xiàn),當階數(shù)不同時,取不同的步長也可以達到較好的濾波效果.當階數(shù)增大時,為了得到較好的濾波效果,步長也將隨之增大.通過仿真我們發(fā)現(xiàn),自適應濾波器在保證快速濾波的前提下,能有效還原信號的輪廓,將誤差降到最小.

4 觀測數(shù)據(jù)驗證

為了檢驗該算法的有效性,本文使用新疆天文臺南山觀測基地觀測的脈沖星數(shù)據(jù)進行了濾波實驗.由于南山26 m射電望遠鏡是單天線觀測,無法獲取實驗中所需的實時參考信號,本文將雙極化信號中的一路作為參考信號,另一路作為主望遠鏡信號,將兩路信號進行對消,驗證自適應濾波器的有效性.實驗采用L波段接收機觀測的脈沖星J0332+5434數(shù)據(jù),數(shù)據(jù)記錄為天文VDIF(VLBI(Very Long Baseline Interferometry)Data Interchange Format)格式,信號帶寬為256 MHz,中心頻率為1428 MHz.數(shù)據(jù)處理時首先將記錄的兩個文 件snap2input02020-11-0710-46-36(極 化1)和snap2input12020-11-0710-46-34(極化2)分別作為主通道和參考通道的信號輸入自適應濾波器,選擇優(yōu)化后的濾波參數(shù)對其進行濾波,階數(shù)為2,步長為0.00015.然后利用VLBI處理軟件DIFX(Distributed FX Correlator)對數(shù)據(jù)進行折疊,處理結果如圖8所示.從上至下,第1、2幅圖分別為未經(jīng)濾波的極化1和極化2脈沖星輪廓,第3幅為采用自適應濾波器處理后的輪廓,具有相關性的兩個輪廓基本被消除.由此可推斷,如將該方法應用至實際觀測中,參考通道與主通道中具有相關輪廓的信號均會被消除,即干擾信號被消除,因為干擾信號將同時出現(xiàn)在參考通道和主通道中.而天文信號將被保留,因為微弱的天文信號僅存在于主通道中.

本文算法還使用Parkes 64 m射電望遠鏡超寬帶接收機經(jīng)系統(tǒng)RFI處理后的測試數(shù)據(jù)進行了算法性能驗證.射電望遠鏡觀測脈沖星J0332+5434,3 m口徑的小天線指向固定干擾源(4G塔),同時記錄數(shù)據(jù).觀測系統(tǒng)將主信道和參考信道的數(shù)據(jù)交替記錄在一個PSRDada格式的文件中,主信道和參考信道的數(shù)據(jù)都是16比特復數(shù)信號(實部和虛部各16比特),帶寬為128 MHz.

首先,利用天文信號處理軟件(DSPSR、PSRCHIVE等)進行消色散及折疊,判斷記錄的數(shù)據(jù)中是否存在脈沖星信號.折疊之后,在主信道中可以看到脈沖星J0332+5434的清晰輪廓,由于參考天線增益低、口徑小、靈敏度低,在實際觀測中無法收到微弱的天文信號,因此在參考信道中看不到脈沖星信號輪廓.然后,使用自適應濾波器進行RFI消除測試,數(shù)據(jù)處理時先將.dada文件中的主信道和參考信道分離出來,再將信號輸入至本模型進行濾波處理,參數(shù)設置為階數(shù)6、步長0.00006.圖9顯示了濾波前后的信號頻域圖,從上至下依次為:主天線信號(Parkes 64 m)、參考天線信號和濾波后的信號.從圖中可以看出,位于0點位置的強干擾信號在經(jīng)過本濾波模型后得到了明顯的削弱.

5 討論與總結

天文信號是一種極其微弱的天體輻射,為了觀測到清晰可見的天文信號需要使用高靈敏度的觀測設備.天線接收面積對觀測靈敏度的影響很大,因此一般使用大口徑射電望遠鏡進行天文觀測.

圖8 濾波前后的輪廓圖Fig.8 The profile figure before and after filtered

圖9 濾波前后的頻域圖Fig.9 Frequency domain figure before and after filtered

RFI主要來自于人類的通信活動和日常生活,例如無線電廣播、電視信號、3G和4G通信信號等,這些干擾信號都很強,通過小天線接收.

主天線和參考天線都接收到天文信號的概率較低,主要原因如下:

(1)主天線與參考天線的指向不同.在觀測過程中,射電望遠鏡跟蹤天文信號源,而參考天線固定指向RFI源,二者指向重合的概率比較低;

(2)參考天線一般采用低增益小口徑天線,靈敏度低,而絕大多數(shù)天文信號都非常微弱,能接收到天文信號的概率非常低.

以Parkes 64 m射電望遠鏡為例,在測量中采用了3 m口徑的參考天線和常溫接收機,信號帶寬為128 MHz,系統(tǒng)噪聲溫度為120 K.對于目前天文上探測到的信號強度極高的快速射電暴(FRB)來說,參考天線的探測靈敏度為:

其中,Smin是在給定閾值γ0下的可探測流量密度,β是數(shù)字化因子,Tsys是系統(tǒng)溫度,G是望遠鏡增益,BW是帶寬,NP是偏振數(shù)[14].當天線效率取60%時,GDPFU(Degree Per Flux Unit)計算后為0.0015 K·Jy?1,τ為FRB信號的爆發(fā)持續(xù)時間,取3 ms計算,探測閾值取7時的最小可探測流量密度Smin約為667 Jy,通過查詢FRB目錄網(wǎng)站1https://frbcat.org/,目前探測到的FRB很少有超過這個數(shù)值的.此外,FRB信號出現(xiàn)在特定天區(qū)的特定時刻,而參考天線固定指向某一特定方向,因此,參考天線接收到天文信號的概率非常低.

本文分析了自適應濾波器的結構和算法,根據(jù)自適應干擾對消原理,設計并建立了適用于射電天文的自適應濾波器仿真模型,通過大量仿真實驗得出了算法參數(shù)與濾波器性能之間的關系,獲取了較優(yōu)的參數(shù)模型,并將南山26 m射電望遠鏡和Parkes 64 m射電望遠鏡的觀測數(shù)據(jù)輸入至該模型,驗證了該自適應濾波器的有效性.自適應濾波器RFI消除方法能濾除干擾信號而保留天文信號,是其他的RFI消除方法所不具備的.由于缺乏參考天線以及相關信號采集設備,本文僅對該方法進行了仿真和部分驗證,作為啟發(fā)性工作,未來將持續(xù)開展相關研究,以期將該方法應用至南山26 m射電望遠鏡、即將建設的新疆110 m大口徑全可動射電望遠鏡以及國內(nèi)其他大型射電望遠鏡.