劉嘯星，蔡體菁，金偉明

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

0 引言

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

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

1 FIR低通濾波器

FIR濾波器的系統(tǒng)函數(shù)：

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

(1)

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

(2)

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

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

2 濾波參數(shù)確定和窗函數(shù)選取

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

(3)

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

(4)

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

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

(5)

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

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

(6)

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

3) 海明窗又稱為改進(jìn)的升余弦窗。當(dāng)將漢寧窗改進(jìn)為海明窗時，旁瓣可更小，長度為N的海明窗函數(shù)：

(7)

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

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

(8)

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

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

表1 窗函數(shù)的基本參數(shù)

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

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

3 試驗結(jié)果比較

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

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

圖1 第一組數(shù)據(jù)濾波結(jié)果

圖2 第二組數(shù)據(jù)濾波結(jié)果

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

圖3 第一組數(shù)據(jù)測區(qū)截圖

圖4 第二組數(shù)據(jù)測區(qū)截圖

對以上兩組測區(qū)數(shù)據(jù)交叉點不符值進(jìn)行計算，表2、3為4種不同濾波窗函數(shù)交叉點不符值的統(tǒng)計結(jié)果。根據(jù)交叉點不符值來選取合適的濾波窗函數(shù)。

表2 第一組數(shù)據(jù)交叉點不符值統(tǒng)計結(jié)果 單位：mGal

表3 第二組數(shù)據(jù)交叉點不符值統(tǒng)計結(jié)果 單位：mGal

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

在第一組數(shù)據(jù)的測區(qū)內(nèi)選取一組南北方向的重復(fù)測線(見圖5)，計算各窗函數(shù)濾波結(jié)果的重復(fù)線內(nèi)符合精度，其統(tǒng)計結(jié)果如表4所示，漢寧窗、海明窗和布萊克曼窗仍表現(xiàn)一致，均符合海洋重力測量數(shù)據(jù)處理的精度需求，其中漢寧窗略優(yōu)。

圖5 第一組測區(qū)重復(fù)測線示意圖

表4 重復(fù)線內(nèi)符合精度統(tǒng)計結(jié)果

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

4 結(jié)束語

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