韋建成，肖 云，王 利，孟 寧，鄒嘉盛

1.長安大學(xué)地質(zhì)工程與測繪學(xué)院，陜西 西安，710054;

2.地理信息工程國家重點實驗室，陜西 西安，710054;

3.西安測繪研究所，陜西 西安，710054

1 引 言

海洋重力測量是在運(yùn)動狀態(tài)下實施重力測量，受到海浪起伏、航速變化、儀器振動以及海風(fēng)和海流等因素的影響，觀測值中不可避免地含有很多測量誤差。重力觀測值會受到水平加速度、垂直加速度、測量船轉(zhuǎn)彎、交叉耦合以及E?tv?s效應(yīng)等多項干擾加速度的影響，一些干擾項的量級往往比實際重力加速度的量級大得多，為了得到高精度重力測量結(jié)果，必須從重力觀測值中剔除這些干擾項的影響，最大程度地還原真實的重力異常信號。

對于消除重力測量誤差，可以從硬件和軟件兩方面入手。在硬件方面，通過增加合適的阻尼系統(tǒng)，消除垂直方向上干擾加速度對重力觀測值的影響。在軟件方面，可以通過設(shè)計合理的數(shù)字濾波器，對重力測量信號進(jìn)行濾波，剔除高頻噪聲，最大程度上還原真實重力信號。相比于增加儀器阻尼的方法，采用數(shù)字濾波器的方法更易實現(xiàn)。在硬件改進(jìn)無法滿足要求時，可以采用濾波的方法來提高測量精度。目前應(yīng)用于海洋重力測量數(shù)據(jù)處理的濾波器種類繁多，主要有RC(Resistor-Capacitor)濾波器、FIR(Finite Impulse Response)濾波器和IIR(Infinite Impulse Response)濾波器等。RC濾波器存在相位延遲，導(dǎo)致其靈敏度較低;IIR濾波器的典型代表巴特沃斯(Butterworth)低通濾波器也在一些海洋重力測量實踐中得到成功應(yīng)用，但這種濾波器不具有線性相位，不能保證濾波前后波形的一致性，并且還存在潛在的不穩(wěn)定性;FIR濾波器的優(yōu)點是具有精確的線性相位而且系統(tǒng)總是穩(wěn)定的。海洋重力測量中采集數(shù)據(jù)的傳感器眾多，數(shù)據(jù)處理中就需要保持濾波前后數(shù)據(jù)的形狀相一致，這就要求所設(shè)計的濾波器具有精確線性相位，F(xiàn)IR濾波器特性滿足此要求。第二代L&R?？罩亓x已經(jīng)開始使用Blackman濾波器(屬于FIR濾波器)替代傳統(tǒng)的RC濾波器。此外，測量數(shù)據(jù)中不可避免地含有粗差，若粗差未能完全剔除，則遞歸運(yùn)算中的粗差累積影響勢必比非遞歸運(yùn)算大，因此，F(xiàn)IR濾波器輸出結(jié)果較IIR濾波器更加穩(wěn)定[1～3]。窗函數(shù)法設(shè)計的濾波器應(yīng)用于海洋重力數(shù)據(jù)處理的文獻(xiàn)相對較多，并且這種方法在實際應(yīng)用中也取得了很好的效果，詳見文獻(xiàn)[1]?；谧罴岩恢卤平ㄔO(shè)計的濾波器在航空重力測量中得到了成功的應(yīng)用，但在海洋重力測量數(shù)據(jù)處理方面還未見有相關(guān)文獻(xiàn);頻率采樣法設(shè)計的濾波器在海洋重力測量數(shù)據(jù)處理中也少有應(yīng)用。所以，本文著重討論后兩種FIR數(shù)字濾波器在海洋重力測量數(shù)據(jù)處理中的設(shè)計與應(yīng)用。

本文基于最佳一致逼近法和頻率采樣法設(shè)計了兩種FIR數(shù)字低通濾波器，并以窗函數(shù)濾波器結(jié)果為依據(jù)，采用模擬數(shù)據(jù)來驗證濾波器的性能。對海洋重力實測數(shù)據(jù)進(jìn)行了濾波處理，詳細(xì)對比了兩種濾波器的效果。

2 數(shù)據(jù)模型

2.1 海洋標(biāo)量重力測量

海洋重力測量的基本模型為：

式中，δg為測點重力擾動;gb代表碼頭處的重力值，由陸地重力儀從重力基點聯(lián)測得到;fZ、fZ0為比力觀測值及其初值;aU為載體垂直加速度改正;δaE為厄特弗斯改正;δaH為水平加速度改正;δaF為空間改正;δaK為零點漂移改正;γ0為橢球面上的正常重力。以上各項改正中均含有高頻誤差，需要進(jìn)行濾波處理。

2.2 FIR低通濾波算法

FIR濾波可以表示為一個差分方程：

式中，x(n)是輸入信號;y(n)是濾波后輸出信號;h(k)是濾波器系數(shù);N是濾波器的長度。因而，濾波的關(guān)鍵在于如何求解出合適的濾波器系數(shù)。

2.2.1 窗函數(shù)法

窗函數(shù)法的設(shè)計思想是用一個有限時寬的h(n)逼近理想的低通濾波器的脈沖響應(yīng)hd(n)，使其頻率響應(yīng)H(ejω)逼近理想低通濾波器的頻率響應(yīng)Hd(ejω)。窗函數(shù)法設(shè)計過程中，h(n)和hd(n)的關(guān)系可以表示為：

式中，w(n)為窗函數(shù)。由于Hanning窗具有較快的衰減速度，選擇Hanning窗來設(shè)計濾波器，實際編程中令主瓣寬度的一半等于實際的歸一化截止頻率來估算濾波器長度。Hanning窗的具體參數(shù)可查閱相關(guān)信號處理書籍。

2.2.2 最佳一致逼近法

最佳一致逼近法也稱為Chebyshev(切比雪夫)逼近法，是在所關(guān)注的頻段內(nèi)使誤差函數(shù)

較為均勻一致，通過選擇合適的H(ejω)和權(quán)函數(shù)W(ejω)，使 E(ejω) 的值達(dá)到最小。權(quán)函數(shù)可通過下式進(jìn)行構(gòu)造

式中，ωP是通帶截止頻率;ωs是阻帶截止頻率;δ1是通帶波紋峰值;δ2為阻帶波紋峰值，它們可通過通帶最大衰減Ap、阻帶最小衰減As換算得到，具體換算關(guān)系可參考文獻(xiàn)[4]。實際編程中可利用Remez交換算法迭代求解出濾波器系數(shù)。

2.2.3 頻率采樣法

頻率采樣法是從頻域出發(fā)，在頻域直接設(shè)計。首先根據(jù)頻域的采樣定理，以等頻率間隔對Hd(ejω) 取樣得到 H(k)：

式中，H(k)為理想低通濾波器的沖擊響應(yīng)，其他符號意義同前。

然后由H(k)通過離散傅里葉逆變換(Inverse Discrete Fourier Transform，IDFT)得到濾波器序列h(n)：

3 算例分析

3.1 垂直加速度模擬及濾波器參數(shù)設(shè)置

對于垂直干擾加速度，由波浪運(yùn)動所引起的測量船垂直方向加速度的周期一般為5～10s，且以6s居多，相對于變化緩慢的重力值來說，該干擾信號屬于高頻擾動。船舶運(yùn)動造成的能量大多集中于1/20～1/3Hz的頻率范圍內(nèi)，而重力異常的頻率一般小于1/120Hz。海況平靜時海浪產(chǎn)生的垂直加速度在15Gal，惡劣時達(dá)到300Gal[5]。

假設(shè)海浪引起的垂直加速度信號的幅值為10mGal，頻率分布在處;垂直干擾加速度信號由3個高頻信號疊加合成，其幅值為50Gal，頻率為z，則量測加速度信號為：(采樣頻率Fs為1Hz)

垂直擾動加速度信號為：

為了檢測濾波器在消除垂向高頻擾動加速度的效果，通過MATLAB軟件進(jìn)行仿真。濾波過程中還存在相位延遲，為此文中采取先將數(shù)據(jù)按順序濾波，然后將所得結(jié)果逆轉(zhuǎn)后再反向濾波的策略，這樣就可實現(xiàn)零相位濾波。為便于后文陳述，這里把采用Hanning窗、最佳一致逼近法和頻率采樣法設(shè)計的FIR低通數(shù)字濾波器分別稱為濾波器1、濾波器2和濾波器3，圖1給出了這三種濾波器的實際脈沖響應(yīng)和幅頻特性曲線，三種濾波器對應(yīng)的濾波器長度N分別為267、221、259階。從圖中可直觀看出，三種濾波器在通帶具有較好的低通性能，在阻帶具有較大的衰減，能夠抑制窗口以外的信號。

3.2 濾波器性能檢測

為了對比濾波器對弱隨機(jī)噪聲和高頻擾動噪聲的濾波效果，本文采取兩種方案對3.1節(jié)模擬的加速度信號進(jìn)行濾波處理。方案一：在量測信號的基礎(chǔ)上添加0.5倍幅值大小的隨機(jī)噪聲，然后對信號進(jìn)行濾波處理;方案二：在量測信號的基礎(chǔ)上添加高頻擾動噪聲，然后對信號進(jìn)行濾波處理。兩種方案的信號時域波形如圖2所示。

圖1 三種FIR低通濾波器的脈沖響應(yīng)和幅頻響應(yīng)

圖2 兩種方案的信號時域波形

由于濾波過程還存在邊界效應(yīng)，需要對數(shù)據(jù)進(jìn)行截斷。本文以濾波器長度N為截斷長度在濾波后數(shù)據(jù)的端部進(jìn)行截斷處理，其起始端點比量測加速度序列的端點落后N個歷元，同時相應(yīng)的序列長度比量測加速度序列短2N個歷元[6]。

為了比較濾波后加速度信號與量測加速度信號的差異，以量測加速度信號為參考值，在濾波值與參考值之間求差。以窗函數(shù)濾波器結(jié)果為依據(jù)，對采用最佳一致逼近法和頻率采樣法設(shè)計的濾波器進(jìn)行對比。對圖2中信號的精度進(jìn)行了統(tǒng)計，其結(jié)果見表1。方案一所得結(jié)果見圖3和表2。方案一結(jié)果顯示，采用Hanning窗、最佳一致逼近法和頻率采樣法設(shè)計的低通濾波器所得加速度序列精度分別為 0.777mGal、0.841mGal和0.914mGal，濾波器2的效果略好于濾波器3，這說明本文設(shè)計的兩種低通濾波器能夠從含有弱隨機(jī)噪聲的加速度信號中提取出原始加速度信號。

表1 兩種方案加速度統(tǒng)計結(jié)果

表2 方案一濾波后加速度與參考值之差的相關(guān)統(tǒng)計信息

圖3 方案一信號截斷后所得濾波值與參考值的差

圖4描述了高頻擾動噪聲信號經(jīng)三種濾波器濾波后所得結(jié)果相對于參考值的變化情況，表3給出了相關(guān)統(tǒng)計信息。由表3可知，濾波器2的結(jié)果也稍好于濾波器3，這與表2得出結(jié)論一致。由于基于最佳一致逼近法比采用頻率采樣法設(shè)計濾波器更容易實現(xiàn)，并且其精度較高，因而最佳一致逼近法設(shè)計的低通濾波器在航空重力測量中得到了廣泛的應(yīng)用。濾波器1在方案一和方案二中都表現(xiàn)出了最好的結(jié)果，源于其較快的衰減速度。這說明本文設(shè)計的低通數(shù)字濾波器，質(zhì)量令人滿意能夠從含有高頻擾動噪聲的信號中還原真實信號。

表3 方案二濾波后加速度與參考值之差的相關(guān)統(tǒng)計信息

為進(jìn)一步分析在實際海洋重力測量數(shù)據(jù)中的應(yīng)用效果，采集了捷聯(lián)式船載重力儀在停泊狀態(tài)的下的重力觀測數(shù)據(jù)，假設(shè)重力異常值不改變，則濾波后重力異常殘差為實際的重力異常誤差。圖5給出了三種濾波器的濾波結(jié)果，為了對比的公正性，對相同時間段的濾波結(jié)果進(jìn)行了統(tǒng)計，結(jié)果見表4。

圖4 方案二信號截斷后所得濾波值與參考值的差

圖5 濾波后的重力異常

表4 濾波后重力異常與參考值之差的相關(guān)統(tǒng)計信息

從表4可看出，濾波器1、2、3的濾波結(jié)果的標(biāo)準(zhǔn)差分別為 0.388mGal、0.386mGal、0.389mGal，標(biāo)準(zhǔn)差都在1mGal以內(nèi)。濾波器2和濾波器3的結(jié)果與濾波器1結(jié)果吻合較好，這說明本文設(shè)計的兩種低通濾波器可以將高頻擾動噪聲分離出來，對濾波后的時間序列進(jìn)行邊界效應(yīng)改正后可以提取出高精度的重力異常信息。

4 結(jié) 語

海洋重力測量數(shù)據(jù)處理中，由于高頻擾動噪聲的存在，需要從重力觀測數(shù)據(jù)中還原出真實重力值，低通濾波就顯得尤為重要?；诖?，本文設(shè)計了兩種適用于海洋重力測量的FIR低通數(shù)字濾波器。以普遍應(yīng)用的窗函數(shù)法設(shè)計的濾波器為依據(jù)，通過仿真分析，驗證了最佳一致逼近法和頻率采樣法設(shè)計的濾波器的性能，結(jié)果表明本文設(shè)計的兩種濾波器具有較好的低通濾波性能。對船載實測重力數(shù)據(jù)濾波后，兩種濾波器都能達(dá)到和窗函數(shù)濾波器相近的精度;同時發(fā)現(xiàn)在相同設(shè)計指標(biāo)下，基于最佳一致逼近法設(shè)計的濾波器低通濾波性能略優(yōu)于采用頻率采樣法設(shè)計的濾波器，而且在達(dá)到相近的濾波效果時，基于最佳一致逼近法設(shè)計的濾波器所需的長度要小于采用頻率采樣法設(shè)計的濾波器的長度，在數(shù)據(jù)截短時會保留更多的有效觀測信息。當(dāng)測線上數(shù)據(jù)量不是很多時，采用最佳一致逼近法設(shè)計的濾波器要比基于頻率采樣法設(shè)計的濾波器更優(yōu)。

[1]歐陽永忠.?？罩亓y量數(shù)據(jù)處理關(guān)鍵技術(shù)研究[J].測繪學(xué)報，2014，43(4)：435-435.

[2]孫中苗，夏哲仁.FIR低通差分器的設(shè)計及其在航空重力測量中的應(yīng)用[J].地球物理學(xué)報，2000，43(6)：850-855.

[3]王靜波.航空重力測量數(shù)據(jù)處理方法技術(shù)研究[D].北京：中國地質(zhì)大學(xué)，2010.

[4]樓順天，李博菡.基于 MATLAB的系統(tǒng)分析與設(shè)計——信號處理[M].西安：西安電子科技大學(xué)出版社，1998.

[5]劉鳳鳴.海洋重力測量數(shù)據(jù)實時處理技術(shù)研究[D].哈爾濱：哈爾濱工程大學(xué)，2008.

[6]周波陽，羅志才，寧津生等.航空矢量重力測量中有限沖激響應(yīng)低通數(shù)字濾波器的設(shè)計[J].武漢大學(xué)學(xué)報·信息科學(xué)版，2015，40(6)：772-778.