基于自適應(yīng)濾波的干擾消除方法研究?

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

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

(2中國(guó)科學(xué)院新疆天文臺(tái) 烏魯木齊 830011)

(3中國(guó)科學(xué)院大學(xué) 北京 100049)

(4新疆微波技術(shù)重點(diǎn)實(shí)驗(yàn)室 烏魯木齊 830011)

1 引言

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

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

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

2 研究?jī)?nèi)容

2.1 自適應(yīng)干擾對(duì)消模型

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

圖1 自適應(yīng)干擾對(duì)消原理圖Fig.1 Schematic diagram of adaptive interference cancellation

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

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

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

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

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

2.2 自適應(yīng)濾波器結(jié)構(gòu)

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

濾波器的輸出為:

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

圖2 FIR濾波器結(jié)構(gòu)Fig.2 The structure of FIR filter

2.3 自適應(yīng)濾波器算法

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

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

最陡下降法迭代計(jì)算權(quán)矢量的公式為:

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

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

算法步驟:

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

(2)對(duì)時(shí)間序列每一時(shí)刻進(jìn)行以下計(jì)算:

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

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

3 仿真與分析

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

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

圖3 輸入信號(hào)Fig.3 The input signal

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

當(dāng)階數(shù)一定,改變步長(zhǎng)時(shí),由圖4可知,計(jì)算結(jié)果隨步長(zhǎng)增加收斂速度變快,步長(zhǎng)越短的計(jì)算時(shí)間越長(zhǎng).濾波性能隨參數(shù)的改變?cè)趫D5中可以看到,誤差的收斂速度與濾波信號(hào)的收斂速度同步,但濾波效果隨步長(zhǎng)增大而變差.

圖5 2階時(shí)不同步長(zhǎng)的誤差(上: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 步長(zhǎng)為0.0001時(shí)不同階數(shù)的濾波信號(hào)(上: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 步長(zhǎng)為0.0001時(shí)不同階數(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)

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

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

4 觀測(cè)數(shù)據(jù)驗(yàn)證

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

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

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

5 討論與總結(jié)

天文信號(hào)是一種極其微弱的天體輻射,為了觀測(cè)到清晰可見(jiàn)的天文信號(hào)需要使用高靈敏度的觀測(cè)設(shè)備.天線接收面積對(duì)觀測(cè)靈敏度的影響很大,因此一般使用大口徑射電望遠(yuǎn)鏡進(jìn)行天文觀測(cè).

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

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

RFI主要來(lái)自于人類(lèi)的通信活動(dòng)和日常生活,例如無(wú)線電廣播、電視信號(hào)、3G和4G通信信號(hào)等,這些干擾信號(hào)都很強(qiáng),通過(guò)小天線接收.

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

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

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

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

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

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