海洋重力數據處理中的正反FIR濾波算法

蔡體菁,邵錦江,胡嘯林

(東南大學 微慣性儀表與先進導航技術教育部重點實驗室,江蘇 南京 210096)

0 引言

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

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

1 線性相位FIR濾波器設計

1.1 基本原理與結構

FIR濾波器的卷積型表達式為

(1)

式中:N為FIR濾波器單位沖激響應h(k)的有限長度;階數M=N-1。其系統函數H(z)在|z|>0處收斂,z平面內只有零點,N-1重極點在z=0處,因此是一個因果穩定的系統;無輸出到輸入的反饋,一般為非遞歸結構。

文獻[5]指出單位沖激響應序列是偶對稱或奇對稱(即h(n)=±h(N-1-n))的濾波器具有線性相位,反之亦然。由此設計的線性相位網絡結構濾波器,最多只需要進行N/2次乘法,降低了濾波器的計算復雜度。當N為奇數,h(n)為偶對稱時,h(n)的相位函數和群延遲為

(2)

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

1.2 基于窗函數的FIR濾波器設計方法

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

(3)

對應的單位沖激響應[7-8]為

(4)

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

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

(5)

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

(6)

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

2 實時正反FIR濾波算法設計

2.1 實時正反FIR濾波公式推導

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

(m≥0)

(7)

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

(8)

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

(9)

式中m不小于j,根據多重求和的運算性質,有:

(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)的運算,由:

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

(16)

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

(17)

計算可知:

(18)

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

(19)

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

(20)

(21)

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

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

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

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

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

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

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

2.3 正反FIR實時濾波方案設計

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

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

2) 調參得到正反FIR濾波器的沖激響應h′1(n),將h′1(0),h′1(1),…,h′1(N1)存儲到嵌入式計算的內存中(設為數組Resp),方便濾波時調用。

3) 開辟一個長度N1的內存空間(設為循環數組DinW),為當數據傳輸到FIR濾波器,按照滑窗的形式依次進入數組DinW,輸入的第(k+1)個數據覆蓋存儲在DinW[k/N1]處。

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

DinW[(k+i)/N1]

(22)

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

圖3 實時正反FIR濾波算法示意圖

3 船載海洋重力試驗

3.1 海洋重力數據特征

截止頻率影響了濾波后重力異常的分辨率及數據的內符合精度。截止頻率的選擇應該兼顧載體航行速度和重力測量的分辨率,三者之間的關系為

λ=v/(2fc)

(23)

(24)

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

3.2 海洋重力數據處理船載試驗

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

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

用以上濾波器對某次海上實驗采集到重力異常數據進行濾波處理,實時單正向FIR濾波結果如圖4所示,實時和事后正反FIR濾波結果去程重復線如圖5所示,返程重復線如圖6所示。

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

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

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

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

4 結束語

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