海洋重力測量中濾波方法比較研究

劉嘯星，蔡體菁，金偉明

(東南大學 儀器科學與工程學院,江蘇 南京 210096)

0 引言

海洋重力輔助慣性導航系統是工作在強噪聲環境下的動態測量平臺，采集的數據中夾雜了不同頻譜特性的噪聲[1-2]。海洋重力測量的數據中噪聲信號常比重力異常信號高出幾百倍甚至上千倍。此外，重力異常信號主要集中在低頻段，此頻段內的噪聲信號幅值通常低于重力異常幅值，而噪聲信號主要集中在高頻段并占取了大部分的帶寬[3-4]。為了獲得實際的重力異常信號，需要對測量數據進行低通濾波處理。

目前常用的濾波方法有有限脈沖響應(FIR)低通濾波、無限脈沖響應(IIR)低通濾波、Kalman濾波和小波濾波等[5-6]。周薇等[7]針對三軸慣性平臺海空重力儀的特點，應用平滑卡爾曼濾波法對海洋重力測量數據進行了濾波處理，根據實測數據分析，其重復線內符合精度為0.2～0.4 mGal(1 Gal=10-2m/s2),空間分辨率達650 m。郎駿健等[8]基于FIR低通濾波器存在邊緣效應而浪費部分有效數據的問題，提出了一種可抑制邊緣效應的傅里葉基追蹤低通濾波方法，該方法采用凸優化中的內點算法，可將低頻與高頻分量有效分離。根據仿真實驗分析和實測數據驗證，該方法東西測線和南北測線的均方根誤差為0.7 mGal和1.4 mGal。本文研究比較了不同窗函數下FIR濾波器處理海洋重力測量數據的效果。

1 FIR低通濾波器

FIR濾波器的系統函數：

H(z)=b0+b1z-1+…+bNz-N=

(1)

式中：N為FIR濾波器單位脈沖響應h(n)的長度。H(z)在Z平面上有N-1個零點，在原點處有1個N-1重極點，所以H(z)永遠穩定。假設理想低通數字濾波器在通帶內的幅頻特性|Hd(ejω)|=1，相頻特性φ(ω)=0，截止頻率ωc=2πfc，具有線性相位，群時延為α，則對應的時域濾波函數為hd(n)，頻率響應為

(2)

FIR數字濾波器要求用有限長的單位脈沖響應h(n)去逼近無限長的理想濾波器的單位脈沖響應hd(n)。目前常用的方法是使用一個有限長的窗函數序列w(n)來截取hd(n)的主要部分h(n)=hd(n)·w(n),n=0,1,2,…,N-1,通過這種方式得到的頻率響應e-jω近似于理想濾波器頻率響應Hd(e-jω)。線性相位濾波器要求h(n)必須是偶對稱的，對稱中心是其長度的1/2，即h(n)=h(N-1-n)，而且α=(N-1)/2，所以要求窗函數w(n)也必須是中心偶對稱的，即w(n)=w(N-1-n)。

由式(1)、(2)可知，基于窗函數的FIR低通濾波器的設計中，濾波器的截止頻率、濾波長度以及窗函數的選取是關鍵因素。

2 濾波參數確定和窗函數選取

低通濾波器截止頻率fc的選擇應兼顧載體航行速度v和重力測量的分辨率λ，其關系為

(3)

fc一般用濾波周期Tc的倒數表示。根據采樣頻率fT可確定歸一化截止頻率為

(4)

考慮到雖然偶數階濾波器的幅頻響應比奇數階濾波器的幅頻響應好，但會產生非整數的時間延遲，濾波后需要進一步進行內插處理，通常濾波器的長度選擇為奇數。令歸一化截止頻率等于窗函數主瓣寬度的一半，則可估算得到濾波器的長度，但該值只是一個估計值，需要以該值為基礎，根據實際情況進行對比修正。對4種窗函數進行選取。

1) 三角形窗。長度為N的三角窗函數：

(5)

三角形窗的主瓣寬度為8π/N，雖然三角形窗的主瓣寬度大于矩形窗的主瓣寬度(4π/N)，但三角形窗的頻譜密度函數永遠是正值。

2) 漢寧窗又稱為升余弦窗。長度為N的漢寧窗函數：

(6)

以上兩部分之和使邊瓣互相抵消，能量更集中在主瓣，邊瓣衰減速度較快。但漢寧窗的主瓣寬度比矩形窗的主瓣寬度大1倍，即為8π/N。

3) 海明窗又稱為改進的升余弦窗。當將漢寧窗改進為海明窗時，旁瓣可更小，長度為N的海明窗函數：

(7)

海明窗可將99.963%的能量集中在窗譜的主瓣內。與漢寧窗相比，主瓣寬度都是8π/N，邊瓣的最大峰值變小了，小于主瓣峰值的1%，但是邊瓣衰減速度相對較慢。

4) 布萊克曼窗又稱為二階2升余弦窗。長度為N的布萊克曼窗函數：

(8)

布萊克曼窗加上了余弦的二次諧波分量來進一步抑制旁瓣峰值，此時的主瓣寬度為12π/N。

窗函數的性能主要由主瓣寬度、旁瓣最大峰值及旁瓣漸近衰減速度來確定。理想的FIR低通濾波器應具有最小的主瓣寬度、最小的旁瓣峰值和最大的旁瓣衰減速度。在實際工程中不太符合理論要求，因此需要綜合截止頻率、濾波長度和窗函數的參數指標來確定濾波器的設計。4種窗函數的性能指標如表1所示。

表1 窗函數的基本參數

根據對上述4種窗函數的頻率特性和窗譜性能指標的描述，三角形窗、漢寧窗和海明窗具有相同的主瓣寬度，旁瓣峰值依次減小。與三角窗相比，漢寧窗和海明窗的旁瓣衰減速度較大，但三角形窗的旁瓣峰值最大。因此，漢寧窗和海明窗更適合于處理海洋重力測量信號。布萊克曼窗具有最小的旁瓣峰值，旁瓣衰減速度也較快，滿足最重要的處理條件，但主瓣寬度最大。因此,雖然漢寧窗、海明窗和布萊克曼窗均初步滿足海洋重力測量數據的處理要求，但無法僅通過性能指標判斷窗函數的優劣，需要通過實測數據的濾波處理結果進行對比來選取最合適的窗函數。

下面針對相同截止頻率下，基于三角窗、漢寧窗、海明窗和布萊克曼窗4種窗函數設計的濾波器，通過對兩組實測數據的比較與分析，選取最適用于海洋重力測量低通濾波的窗函數。

3 試驗結果比較

選取某兩個海區重力測量數據進行濾波方法比較。載體的航行速度均為6 m/s，采樣頻率fT=1 Hz。根據本次采集數據的空間分辨率的要求及式(3)可得截止頻率為0.003 33 Hz，即濾波周期為300 s，根據式(4)可計算出相應的歸一化截止頻率為0.003 33 Hz。由窗函數的主瓣寬度和歸一化截止頻率可計算出各個濾波器的濾波長度。由于本次采用奇數長度的濾波器，這里三角窗、漢寧窗和海明窗的長度為601，布萊克曼窗的長度為901。分別使用以上各參數對測線數據進行濾波處理，分析不同的窗函數對于海洋重力測量數據處理的適用性，最佳窗函數的選擇根據重復線內符合精度和測線交叉點不符值來進行評估。

圖1、2分別為第一、二組數據濾波結果。由圖1、2可看出，三角窗函數進行濾波后，測線數據仍有較明顯的振蕩，數據中高頻噪聲未能被完全濾除，這驗證了前文所述的三角窗的旁瓣峰值大、旁瓣衰減速度小而不適合處理海洋重力測量數據的結論。漢寧窗、海明窗和布萊克曼窗由于旁瓣峰值小及旁瓣衰減速度快，對測線濾波后的數據相對平滑，可有效地去除高頻噪聲，由圖還可看出三者的濾波效果水平相當。

圖1 第一組數據濾波結果

圖2 第二組數據濾波結果

基于以上4種窗函數濾波后的結果，分別對兩組數據測區的航跡進行了截取，如圖3、4所示。由圖可看出，第一組測區的測線呈均勻的縱橫分布，其中東西測線5條，南北測線6條，共計算交叉點30個。第二組測區的測線分布不均勻且形狀不規則，此處截取東西測線4條，南北測線2條，共計算交叉點不符值8個。

圖3 第一組數據測區截圖

圖4 第二組數據測區截圖

對以上兩組測區數據交叉點不符值進行計算，表2、3為4種不同濾波窗函數交叉點不符值的統計結果。根據交叉點不符值來選取合適的濾波窗函數。

表2 第一組數據交叉點不符值統計結果 單位：mGal

表3 第二組數據交叉點不符值統計結果 單位：mGal

由表2、3可看出，雖然三角窗、漢寧窗和海明窗的濾波長度一致，但三者中三角窗的濾波效果最差，其中誤差分別為0.73 mGal和1.41 mGal。漢寧窗、海明窗和布萊克曼窗的濾波效果相當，其中漢寧窗和海明窗效果最好，其中誤差均為0.47 mGal和0.41 mGal。布萊克曼窗的濾波長度最長，與另外兩個窗函數的濾波結果相比，其在濾波過程中舍棄了測線頭尾各一個濾波周期的數據，造成了有效數據的浪費。

在第一組數據的測區內選取一組南北方向的重復測線(見圖5)，計算各窗函數濾波結果的重復線內符合精度，其統計結果如表4所示，漢寧窗、海明窗和布萊克曼窗仍表現一致，均符合海洋重力測量數據處理的精度需求，其中漢寧窗略優。

圖5 第一組測區重復測線示意圖

表4 重復線內符合精度統計結果

綜上所述，基于漢寧窗、海明窗和布萊克曼窗設計的濾波器對于船測重力數據的濾波效果都適用，綜合考慮有效數據的保留和交叉點不符值的結果，漢寧窗和海明窗對于海洋重力測量數據的處理是一個理想的選擇。

4 結束語

采用低通濾波器來提取海洋重力測量數據中混合在高頻噪聲中的低頻重力異常信息。本文給出了FIR低通數字濾波器截止頻率和濾波器長度的確定方法，結合兩組實測的海洋重力數據，通過對不同窗函數的濾波結果計算其重復線內符合精度和交叉點不符值，漢寧窗、海明窗和布萊克曼窗的計算精度均優于0.5 mGal，符合海洋重力測量數據處理的要求，綜合考慮有效數據的保留和交叉點不符值的結果，選擇漢寧窗和海明窗設計的FIR低通數字濾波器來提取海洋重力場信息。