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

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

2.地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安，710054;

3.西安測(cè)繪研究所，陜西 西安，710054

1 引 言

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

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

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

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

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

海洋重力測(cè)量的基本模型為：

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

2.2 FIR低通濾波算法

FIR濾波可以表示為一個(gè)差分方程：

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

2.2.1 窗函數(shù)法

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

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

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]。實(shí)際編程中可利用Remez交換算法迭代求解出濾波器系數(shù)。

2.2.3 頻率采樣法

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

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

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

3 算例分析

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

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

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

垂直擾動(dòng)加速度信號(hào)為：

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

3.2 濾波器性能檢測(cè)

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

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

圖2 兩種方案的信號(hào)時(shí)域波形

由于濾波過程還存在邊界效應(yīng)，需要對(duì)數(shù)據(jù)進(jìn)行截?cái)?。本文以濾波器長(zhǎng)度N為截?cái)嚅L(zhǎng)度在濾波后數(shù)據(jù)的端部進(jìn)行截?cái)嗵幚?，其起始端點(diǎn)比量測(cè)加速度序列的端點(diǎn)落后N個(gè)歷元，同時(shí)相應(yīng)的序列長(zhǎng)度比量測(cè)加速度序列短2N個(gè)歷元[6]。

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

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

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

圖3 方案一信號(hào)截?cái)嗪笏脼V波值與參考值的差

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

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

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

圖4 方案二信號(hào)截?cái)嗪笏脼V波值與參考值的差

圖5 濾波后的重力異常

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

從表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é)果吻合較好，這說(shuō)明本文設(shè)計(jì)的兩種低通濾波器可以將高頻擾動(dòng)噪聲分離出來(lái)，對(duì)濾波后的時(shí)間序列進(jìn)行邊界效應(yīng)改正后可以提取出高精度的重力異常信息。

4 結(jié) 語(yǔ)

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

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

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

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

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

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

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