蔡體菁,邵錦江,胡嘯林

(東南大學(xué) 微慣性儀表與先進(jìn)導(dǎo)航技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,江蘇 南京 210096)

0 引言

在海洋重力場信號測量時(shí),由于環(huán)境和儀器影響,大量的測量噪聲也會被傳感器所采集。目前在重力數(shù)據(jù)處理上有很多濾波技術(shù),除了基于頻域?yàn)V波的有限沖擊響應(yīng)(FIR)濾波、無限脈沖響應(yīng)(IIR)濾波外,基于特征的小波濾波和基于貝葉斯估計(jì)的Kalman濾波方法[1]。

對重力測量數(shù)據(jù)濾波處理時(shí),由于重力異常信號一般頻率非常低,所以數(shù)字低通濾波器的應(yīng)用比較廣泛。FIR濾波具有線性相位,易于設(shè)計(jì),IIR相位非線性,設(shè)計(jì)較復(fù)雜。因此,常采用FIR低通濾波器去剔除重力測量信號中的高頻噪聲[2-3]。郭圣煥[4]在比較漢寧窗FIR、切比雪夫等紋波FIR和巴特沃斯IIR低通數(shù)字濾波器后發(fā)現(xiàn),得出切比雪夫紋波FIR低通濾波器具有更高的濾波精度。但FIR對于與有用信號頻段接近的噪聲無能為力,并且低通濾波需要一段數(shù)據(jù)才能進(jìn)行一次濾波,因此導(dǎo)致了濾波結(jié)果的滯后性,這是不利于實(shí)際應(yīng)用的。若能實(shí)現(xiàn)全頻段線性相位,則濾波器不會出現(xiàn)非線性失真,這對后續(xù)的濾波組合算法有十分重大的意義。

1 線性相位FIR濾波器設(shè)計(jì)

1.1 基本原理與結(jié)構(gòu)

FIR濾波器的卷積型表達(dá)式為

(1)

式中:N為FIR濾波器單位沖激響應(yīng)h(k)的有限長度;階數(shù)M=N-1。其系統(tǒng)函數(shù)H(z)在|z|>0處收斂,z平面內(nèi)只有零點(diǎn),N-1重極點(diǎn)在z=0處,因此是一個(gè)因果穩(wěn)定的系統(tǒng);無輸出到輸入的反饋,一般為非遞歸結(jié)構(gòu)。

文獻(xiàn)[5]指出單位沖激響應(yīng)序列是偶對稱或奇對稱(即h(n)=±h(N-1-n))的濾波器具有線性相位,反之亦然。由此設(shè)計(jì)的線性相位網(wǎng)絡(luò)結(jié)構(gòu)濾波器,最多只需要進(jìn)行N/2次乘法,降低了濾波器的計(jì)算復(fù)雜度。當(dāng)N為奇數(shù),h(n)為偶對稱時(shí),h(n)的相位函數(shù)和群延遲為

(2)

即y(n)表示(n-α)時(shí)刻數(shù)據(jù)的濾波結(jié)果,經(jīng)分析可知FIR濾波器的相位延遲等于濾波器階數(shù)的1/2。FIR濾波器的相位特性主要受h(n)的奇偶對稱性影響,同時(shí)也受H(ω)影響。如果幅度函數(shù)H(ω)>0,則系統(tǒng)的相位等于-(N-1)/2·ω;如果H(ω)<0,則系統(tǒng)的相位等于-(N-1)/2·ω+π,從這個(gè)意義上來講,這樣的系統(tǒng)并不滿足嚴(yán)格意義上的線性相頻特性[6]。只有當(dāng)對任意ω都有H(ω)>0時(shí),系統(tǒng)才是嚴(yán)格意義上的線性相位系統(tǒng)。

1.2 基于窗函數(shù)的FIR濾波器設(shè)計(jì)方法

給定一個(gè)理想的零相位低通濾波器,設(shè)濾波的截止頻率為ωc,其頻率特性可表示為

(3)

對應(yīng)的單位沖激響應(yīng)[7-8]為

(4)

對應(yīng)FIR濾波器設(shè)計(jì)的目標(biāo)是尋找一個(gè)頻響函數(shù)H(ejω)去逼近Hd(ejω),逼近的方式有頻率采樣法和時(shí)域的窗函數(shù)法。時(shí)域的窗函數(shù)設(shè)計(jì)思路是設(shè)計(jì)有限長的h(n)逼近hd(n)。最簡單的逼近方法是用一個(gè)窗口截取一段hd(n),此時(shí)對應(yīng)的有限長度的FIR沖激響應(yīng)h(n)可以表示為hd(n)和一個(gè)窗函數(shù)w(n)的乘積[9],即:

h(n)=w(n)hd(n)

(5)

窗函數(shù)w(n)除了取最簡單的矩形脈沖函數(shù)外,還可以取其他形式對hd(n)作一定的權(quán)重變化處理,來改善濾波器的特性。此外,H(ejω)的通帶和阻帶范圍內(nèi)的波動幅度取決于窗函數(shù)的傅里葉變換W(ejω)的旁瓣大小和數(shù)量,而過渡帶寬則由其主瓣決定。例如漢寧窗表達(dá)式如下:

(6)

式中RN(n)為矩形窗函數(shù)。經(jīng)計(jì)算可知,漢寧窗阻帶最小衰減為-44 dB,旁瓣峰值為-31 dB,過渡帶寬6.2 π/N。其他的窗函數(shù)還包括布萊克曼窗、三角窗、海明窗等等。實(shí)際應(yīng)用中,根據(jù)允許的過渡帶寬以及阻帶衰減等參數(shù),初步選擇窗函數(shù)的形式以及FIR濾波器的階數(shù),再根據(jù)仿真結(jié)果進(jìn)一步調(diào)整參數(shù)直至滿足設(shè)計(jì)要求。

2 實(shí)時(shí)正反FIR濾波算法設(shè)計(jì)

2.1 實(shí)時(shí)正反FIR濾波公式推導(dǎo)

設(shè)N為奇數(shù),h1(n)為偶對稱序列,用小寫z表示n時(shí)刻輸入數(shù)據(jù)對應(yīng)的濾波輸出結(jié)果。結(jié)合式(1)與群延遲計(jì)算公式可知:

(m≥0)

(7)

z′1[m+(N-1)/2]表示m+(N-1)/2時(shí)刻的數(shù)據(jù)的濾波輸出結(jié)果,當(dāng)濾波輸出數(shù)據(jù)量m≥(N-1)時(shí),對輸入的數(shù)據(jù)用同一濾波器進(jìn)行反向FIR濾波,濾波過程實(shí)際上是一種線性卷積,考慮到相位延遲,設(shè)m時(shí)刻輸入數(shù)據(jù)的濾波結(jié)果為z1(m),由于階數(shù)M=N-1,根據(jù)式(1)形式有:

(8)

將式(7)代入式(8)可得:

(9)

式中m不小于j,根據(jù)多重求和的運(yùn)算性質(zhì),有:

(10)

其中m≥(N-1),再令:

ai,j=h1(j)h1(i)

(11)

xi,j=x(m+i-j)

(12)

A=(ai,j)∈N×N

(13)

X=(xi,j)∈N×N

(14)

則:

z1(m)=Tr(AXT)

(15)

式中Tr()表示求跡,由于h1(n)是偶對稱序列,則ai,j=aM-I,j=ai,M-j=aM-i,M-j=aj,i,為了簡化式(10)的運(yùn)算,由:

xi+δ,j+δ=x(m+i-j)

(16)

可得矩陣X所有行列連續(xù)的子矩陣對角線上的元素均相等,令k=i-j,則k∈{-M≤k≤M,k∈Z},則:

(17)

計(jì)算可知:

(18)

此時(shí)h″1(k)為非因果序列,化簡h″1(k),同時(shí)用(-k+M)替換k可得:

(19)

取n=m+M,取N1=2N-1,有M1=2M,那么:

(20)

(21)

式中M=α1=(N1-1)/2表示濾波器的延遲。至此得到了正反FIR濾波器y1(n)式(1)形式的表達(dá)式,能夠?qū)崟r(shí)濾波,運(yùn)算量從式(10)中的N2次乘法運(yùn)算縮短為2N次乘法運(yùn)算,進(jìn)一步提高了正反FIR濾波的效率。分析沖激響應(yīng)h′1(M1-k)=h′1(k),可知h′1(n)也為偶對稱序列。因此,本節(jié)詳細(xì)推導(dǎo)了一個(gè)窗長N1為奇數(shù),沖激響應(yīng)h′1(n)為偶對稱序列的線性相位實(shí)時(shí)正反FIR濾波器表達(dá)式。

2.2 正反FIR濾波頻率特性分析

采用單正向FIR濾波器y2(n)做為對比,其重啟響應(yīng)h2(n)也為偶對稱序列,且采用與h1(n)相同類型的窗函數(shù)設(shè)計(jì)得到,與h′1(n)長度相同。設(shè)整個(gè)FIR濾波器的截至頻率ωc=0.25π,濾波器長度N1=61,窗函數(shù)類型為漢寧窗,設(shè)計(jì)單正向FIR濾波器y2(n),其幅頻與相頻特性曲線如圖1所示。

圖1 漢寧窗60階ωc=0.25π單正向FIR濾波器特性曲線

設(shè)計(jì)正反FIR濾波器時(shí),考慮到h′1(n)需要與h2(n)序列長度相同,則h1(n)的序列長度N=(N1+1)/2=31,同樣取截止頻率ωc=0.25π,再通過式(18)計(jì)算得到h′1(n),則h′1(n)為長度N1=61的正反FIR濾波器y1(n)沖激響應(yīng)。由于在計(jì)算過程中截止頻率發(fā)生偏移,需要微調(diào)h1(n)對應(yīng)的截止頻率,使得y1(n)的截止頻率與y2(n)的相同,最終得到正反向FIR濾波器y1(n),其幅頻與相頻特性曲線如圖2所示。

圖2 漢寧窗60階ωc=0.25π正反FIR濾波器特性曲線

分析可知,正反FIR濾波最顯著的特點(diǎn)是濾波器在所有的采樣頻譜都滿足嚴(yán)格意義上的線性相頻特性,不會對信號造成非線性失真,有利于后續(xù)對信號的進(jìn)一步處理。其次正反FIR濾波的阻帶衰減大,幾乎是單正向FIR濾波的2倍,需要改進(jìn)的地方是過渡帶寬大,實(shí)際應(yīng)用中要考慮選擇合適的參數(shù)以滿足設(shè)計(jì)要求。

2.3 正反FIR實(shí)時(shí)濾波方案設(shè)計(jì)

實(shí)時(shí)處理算法與后處理算法的區(qū)別:實(shí)時(shí)處理算法輸入數(shù)據(jù)在經(jīng)歷過濾波器的延遲時(shí)間后即刻輸出,而后處理算法則需要讀取全部數(shù)據(jù),處理完再輸出,顯然實(shí)時(shí)處理算法更適合處理船載試驗(yàn)中的海洋重力異常數(shù)據(jù)。

1) 實(shí)時(shí)正反FIR濾波算法首先要根據(jù)濾波器的截止頻率,通帶平穩(wěn)度,允許的過渡帶寬以及阻帶衰減,初步選擇窗函數(shù)與濾波器階數(shù)。

2) 調(diào)參得到正反FIR濾波器的沖激響應(yīng)h′1(n),將h′1(0),h′1(1),…,h′1(N1)存儲到嵌入式計(jì)算的內(nèi)存中(設(shè)為數(shù)組Resp),方便濾波時(shí)調(diào)用。

3) 開辟一個(gè)長度N1的內(nèi)存空間(設(shè)為循環(huán)數(shù)組DinW),為當(dāng)數(shù)據(jù)傳輸?shù)紽IR濾波器,按照滑窗的形式依次進(jìn)入數(shù)組DinW,輸入的第(k+1)個(gè)數(shù)據(jù)覆蓋存儲在DinW[k/N1]處。

4) 將數(shù)組Resp中的數(shù)據(jù)與DinW[(k+1)/N1]～DinW[(k+N1)/N1]組成的新的數(shù)組對應(yīng)位置中的數(shù)據(jù)求和:

DinW[(k+i)/N1]

(22)

即為延遲了(N1-1)/2輸出的濾波結(jié)果。綜上所述,實(shí)時(shí)正反FIR濾波算法的示意圖如圖3所示。

圖3 實(shí)時(shí)正反FIR濾波算法示意圖

3 船載海洋重力試驗(yàn)

3.1 海洋重力數(shù)據(jù)特征

截止頻率影響了濾波后重力異常的分辨率及數(shù)據(jù)的內(nèi)符合精度。截止頻率的選擇應(yīng)該兼顧載體航行速度和重力測量的分辨率,三者之間的關(guān)系為

λ=v/(2fc)

(23)

(24)

式中:v為載體的航行速度;λ為重力測量數(shù)據(jù)的空間分辨率;fc=1/Tc為低通濾波器的截止頻率,Tc為濾波周期。

3.2 海洋重力數(shù)據(jù)處理船載試驗(yàn)

針對相同截止頻率下,基于漢寧窗設(shè)計(jì)的實(shí)時(shí)正反FIR濾波器對實(shí)際數(shù)據(jù)濾波效果進(jìn)行研究分析,其中評價(jià)窗函數(shù)性能通過重復(fù)線內(nèi)符合精度進(jìn)行評估。對比實(shí)驗(yàn)組1為實(shí)時(shí)單正向FIR濾波器,對比實(shí)驗(yàn)組2為離線(事后處理)正反FIR濾波器。

試驗(yàn)選取測線數(shù)據(jù)中載體航行速度為6 m/s,采樣頻率為fT=1 Hz,海洋重力測量要求重力異常數(shù)據(jù)空間分辨率達(dá)0.9 km,根據(jù)式(24)可以得到截止頻率為0.003 33 Hz,即濾波周期為300 s,計(jì)算出相應(yīng)的歸一化截止頻率為0.003 33 Hz。計(jì)算得到FIR濾波器的階數(shù)N1=601,調(diào)參得到漢寧窗單正向FIR濾波器和實(shí)時(shí)在線和事后離線FIR濾波器。

用以上濾波器對某次海上實(shí)驗(yàn)采集到重力異常數(shù)據(jù)進(jìn)行濾波處理,實(shí)時(shí)單正向FIR濾波結(jié)果如圖4所示,實(shí)時(shí)和事后正反FIR濾波結(jié)果去程重復(fù)線如圖5所示,返程重復(fù)線如圖6所示。

圖4 漢寧窗600階fc= 0.003 33 Hz 單正向FIR濾波效果

圖5 漢寧窗在線、離線正反FIR去程重復(fù)線

圖6 漢寧窗在線、離線正反FIR返程重復(fù)線

分析可知單正向FIR濾波在輸入數(shù)據(jù)噪聲方差較大時(shí),效果比較差,而實(shí)時(shí)正反FIR濾波則能很好的抑制輸入噪聲,并且其處理重力異常數(shù)據(jù)的重復(fù)線精度為0.09 mGal,與事后處理基本相當(dāng),符合海洋重力數(shù)據(jù)實(shí)時(shí)處理的要求。

4 結(jié)束語

本文提出一種實(shí)時(shí)正反FIR濾波算法,該算法設(shè)計(jì)出的濾波器在所有的采樣頻譜內(nèi)都滿足嚴(yán)格意義上的線性相頻特性,不會對信號造成非線性失真,有利于后續(xù)對信號的進(jìn)一步處理。該正反FIR濾波器的阻帶衰減大,幾乎是單正向FIR濾波的兩倍,需要改進(jìn)的地方是過渡帶寬大,實(shí)際應(yīng)用中要考慮選擇合適的參數(shù)以滿足設(shè)計(jì)要求。船載試驗(yàn)表明,將該實(shí)時(shí)正反FIR濾波器應(yīng)用到海洋重力數(shù)據(jù)處理領(lǐng)域,濾波效果優(yōu)于傳統(tǒng)的單正向FIR濾波器,與事后正反FIR濾波效果相當(dāng)。