海洋重力測量中低通濾波器的設計與應用

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

1.長安大學地質工程與測繪學院，陜西 西安，710054;

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

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

1 引 言

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

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

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

2 數據模型

2.1 海洋標量重力測量

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

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

2.2 FIR低通濾波算法

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

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

2.2.1 窗函數法

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

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

2.2.2 最佳一致逼近法

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

較為均勻一致，通過選擇合適的H(ejω)和權函數W(ejω)，使 E(ejω) 的值達到最小。權函數可通過下式進行構造

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

2.2.3 頻率采樣法

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

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

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

3 算例分析

3.1 垂直加速度模擬及濾波器參數設置

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

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

垂直擾動加速度信號為：

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

3.2 濾波器性能檢測

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

圖1 三種FIR低通濾波器的脈沖響應和幅頻響應

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

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

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

表1 兩種方案加速度統計結果

表2 方案一濾波后加速度與參考值之差的相關統計信息

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

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

表3 方案二濾波后加速度與參考值之差的相關統計信息

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

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

圖5 濾波后的重力異常

表4 濾波后重力異常與參考值之差的相關統計信息

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

4 結 語

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

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

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

[3]王靜波.航空重力測量數據處理方法技術研究[D].北京：中國地質大學，2010.

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

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

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