鄭 崴，王建玲，陳學(xué)鋒

航空重力垂直加速度改正濾波方法對(duì)比研究

鄭 崴，王建玲，陳學(xué)鋒

（河南工學(xué)院 電子信息工程學(xué)院，河南 新鄉(xiāng) 453003）

根據(jù)航空重力測(cè)量的基本原理，闡述了用于航空重力垂直加速度改正的頻域?yàn)V波器和時(shí)域?yàn)V波器的設(shè)計(jì)方法，以此為依據(jù)設(shè)計(jì)了用于垂直加速度改正的低通濾波器和卡爾曼平滑濾波器，并利用濾波器對(duì)航空重力實(shí)測(cè)數(shù)據(jù)進(jìn)行了處理試驗(yàn)。試驗(yàn)結(jié)果表明：兩種濾波方法都能夠?qū)崿F(xiàn)垂直加速度改正，而時(shí)域?yàn)V波器處理精度較高，與解算參考標(biāo)準(zhǔn)的誤差低于1mGal。

航空重力測(cè)量；垂直加速度改正；低通濾波；卡爾曼濾波平滑

0 引言

航空重力測(cè)量是以飛機(jī)為載體，測(cè)量近地空中重力場(chǎng)信息的動(dòng)態(tài)測(cè)量技術(shù)。它能夠在人力無法達(dá)到的區(qū)域完成高精度、大面積的重力測(cè)量任務(wù)[1-4]。在測(cè)量過程中，它綜合運(yùn)用了航空重力儀、全球定位系統(tǒng)GPS（Global Positioning System）、慣性導(dǎo)航系統(tǒng)INS（Inertial Navigation System）等多種技術(shù)設(shè)備。航空重力儀測(cè)得的數(shù)據(jù)是重力信息與多種擾動(dòng)加速度的總和，為獲得高精度的重力異常信號(hào)，需對(duì)測(cè)量數(shù)據(jù)進(jìn)行相應(yīng)的改正，以消除擾動(dòng)加速度的影響。其中，一部分的擾動(dòng)加速度（如厄特缶斯加速度、水平擾動(dòng)加速度）可以利用較為嚴(yán)密的解析公式通過計(jì)算來消除。而飛機(jī)垂直運(yùn)動(dòng)產(chǎn)生的擾動(dòng)加速度無法直接觀測(cè)，也沒有解析公式可以直接計(jì)算得出，因此，為了實(shí)現(xiàn)垂直加速度改正，需要采用濾波技術(shù)來消除干擾[5-7]。

基于所采用的濾波技術(shù)的差別，垂直加速度改正又可分為基于頻域?yàn)V波器的處理方法和基于時(shí)域卡爾曼濾波器的處理方法。前者是先利用數(shù)字差分對(duì)GPS定位信息進(jìn)行處理以獲得載體的垂直運(yùn)動(dòng)加速度，再將其與航空重力儀測(cè)得的數(shù)據(jù)求差，最后通過低通濾波器獲取高精度的重力異常[8-10]。該種處理方法得以實(shí)現(xiàn)的依據(jù)在于：只有在信號(hào)低頻端很小的范圍內(nèi)，重力信號(hào)的能量會(huì)高于噪聲干擾，因此通過低通濾波器限制信號(hào)帶寬，降低數(shù)據(jù)中的高頻噪聲干擾，可獲取滿足精度要求的重力異常。目前，常用到的處理方法包括無限脈沖響應(yīng)（Infinite Impulse Response, IIR）濾波法、有限沖擊響應(yīng)（Finite Impulse Response, FIR）濾波法、快速傅里葉變換濾波法、小波變換濾波法等[11-13]。

而基于時(shí)域?yàn)V波器的處理方法則指的是基于卡爾曼濾波的處理方法，它是通過利用GPS/INS及穩(wěn)定平臺(tái)的觀測(cè)數(shù)據(jù)對(duì)重力測(cè)量數(shù)據(jù)進(jìn)行修正，在建立系統(tǒng)方程的基礎(chǔ)上，以最小方差準(zhǔn)則結(jié)合觀測(cè)數(shù)據(jù)對(duì)重力異常進(jìn)行估計(jì)，進(jìn)而實(shí)現(xiàn)垂直加速度改正。俄羅斯國(guó)立大學(xué)的Bolotin等人基于航空重力儀沿飛行測(cè)線的動(dòng)力學(xué)方程，建立了基于卡爾曼濾波的垂直加速度改正處理模型，實(shí)驗(yàn)結(jié)果表明該模型可獲得精度在毫伽級(jí)別的重力異常[14-16]。由于種種原因，Bolotin等人只在其相關(guān)文獻(xiàn)資料里介紹了該種方法的設(shè)計(jì)思路，但對(duì)于狀態(tài)空間模型構(gòu)建及相應(yīng)的濾波器設(shè)計(jì)細(xì)節(jié)并未提及。張貴賓教授的課題組通過摸索和鉆研，構(gòu)建了基于重力異常模型的信息提取方法，并設(shè)計(jì)了固定區(qū)間平滑器，進(jìn)一步提高了信息提取的精度[17][18]。

本文在航空重力測(cè)量、重力異常提取基本原理的基礎(chǔ)上，根據(jù)航空重力測(cè)量的基本數(shù)學(xué)模型，依據(jù)不同濾波器處理的特點(diǎn)，分別設(shè)計(jì)了用于垂直加速度改正的頻域?yàn)V波器和時(shí)域卡爾曼平滑濾波器，通過實(shí)測(cè)數(shù)據(jù)驗(yàn)證對(duì)比兩種濾波器所獲取的重力異常的精度，并探討它們的特點(diǎn)。

1 航空重力測(cè)量基本原理

航空重力測(cè)量基本原理是從航空重力儀測(cè)得的比力數(shù)據(jù)（重力加速度與飛機(jī)垂直加速度之和）中減去垂直擾動(dòng)加速度，最終獲取所需的重力異常。依據(jù)牛頓第二定律，航空重力儀傳感器質(zhì)點(diǎn)在垂直方向的方程可表示為[19]：

正常場(chǎng)及高度改正項(xiàng)的計(jì)算則有[19]：

為方便分析，對(duì)式（1）進(jìn)行調(diào)整，得：

2 卡爾曼濾波平滑處理

卡爾曼濾波是指在時(shí)域內(nèi)，依據(jù)構(gòu)造的狀態(tài)空間方程對(duì)所需的狀態(tài)向量進(jìn)行估計(jì)的處理方法。為利用卡爾曼平滑濾波實(shí)現(xiàn)垂直加速度改正，首先需要依據(jù)式（1）建立航空重力場(chǎng)測(cè)量的狀態(tài)空間模型，它由系統(tǒng)方程和量測(cè)方程組成。

在航空重力測(cè)量中，將航空重力儀沿測(cè)線的實(shí)時(shí)高程作為卡爾曼濾波器的量測(cè)方程，有：

在測(cè)量時(shí)，航空重力儀會(huì)因?yàn)榉€(wěn)定平臺(tái)的舒拉或阻尼振蕩而偏離水平方向，可以利用測(cè)得的穩(wěn)定平臺(tái)偏離水平方向的失準(zhǔn)角對(duì)其進(jìn)行修正，則式（4）可寫為：

航空重力測(cè)量中，重力異?？梢钥醋魇茄販y(cè)線分布的隨機(jī)過程，因此用隨機(jī)統(tǒng)計(jì)模型對(duì)其進(jìn)行描述，有[20]：

將式（6）和（7）聯(lián)立，可得：

將式（8）作為卡爾曼濾波的系統(tǒng)方程，將式（5）作為卡爾曼濾波的量測(cè)方程，并用離散系統(tǒng)狀態(tài)空間形式表示為：

3 FIR低通濾波處理

與卡爾曼濾波處理算法不同，基于低通濾波器的垂直加速度改正處理方法是一種較為傳統(tǒng)的處理方法。大量的研究表明，具有嚴(yán)格線性相位特征的FIR濾波器在航空重力數(shù)據(jù)處理中具有較好的效果。在信號(hào)處理技術(shù)高度發(fā)展的今天，濾波器的設(shè)計(jì)已經(jīng)不是復(fù)雜的問題，重要的是如何根據(jù)信號(hào)特性選擇合適的濾波器參數(shù)，因此本部分著重討論載體垂直加速度的確定方法和FIR低通濾波器參數(shù)的選取方法，有關(guān)FIR低通濾波器的設(shè)計(jì)方法可參見相關(guān)文獻(xiàn)[8][22]。

為實(shí)現(xiàn)基于FIR低通濾波的垂直加速度改正，首先要確定載體的垂直運(yùn)動(dòng)速度，慣用方法是對(duì)GPS測(cè)定的速度序列求一次差分運(yùn)算，而GPS測(cè)速方式主要有以下三種[22]：

（1）位置微分法，利用GPS提供的位置數(shù)據(jù)進(jìn)行差分運(yùn)算，從而獲得載體的運(yùn)動(dòng)速度；

（2）利用多普勒原始觀測(cè)值直接計(jì)算載體速度；

（3）利用由載波相位中心差分獲得的多普勒觀測(cè)值導(dǎo)出載體速度。

已有研究指出，三種測(cè)速方法的精度基本相當(dāng)，而位置微分法和載波相位中心差分測(cè)速精度略優(yōu)于多普勒原始測(cè)速的精度[23]。

本文利用GPS提供的垂直速度進(jìn)行一次差分來計(jì)算載體的垂直加速度。由于重力異常信號(hào)處于0.01 Hz的低頻段內(nèi)，因此采用兩點(diǎn)中心差分即可滿足精度要求[22]。兩點(diǎn)中心差分期輸入與輸出關(guān)系有：

其傳遞函數(shù)為：

求得載體垂直加速度后，由于利用式（4）計(jì)算所獲得的結(jié)果中還包含著大量的噪聲干擾，還需利用FIR低通濾波器對(duì)其進(jìn)行處理后，才能夠獲得高精度的重力異常。利用窗函數(shù)法設(shè)計(jì)FIR低通濾波器時(shí)，需要確定三個(gè)參數(shù)：窗函數(shù)的類型、窗函數(shù)的階數(shù)和截止頻率f。具體來說，在選取窗函數(shù)類型時(shí)，需要根據(jù)原始重力異常的頻域特點(diǎn)來選擇滿足要求的窗函數(shù)，在航空重力測(cè)量數(shù)據(jù)處理中，海明窗、漢寧窗、布雷克曼窗和凱澤窗函數(shù)的性能比較優(yōu)越[24]；對(duì)于來說，為保證濾波器擁有較好的幅頻特性，使其能夠準(zhǔn)確地分離噪聲和信號(hào)，應(yīng)選擇較寬的濾波窗口，但這會(huì)導(dǎo)致處理后有效數(shù)據(jù)點(diǎn)的減少，反之如降低濾波窗口寬度可減弱邊界效應(yīng)的影響，但濾波器的幅頻特性便不能保證；對(duì)于截止頻率f的選擇來說，要綜合考慮重力異常精度和分辨率。航空重力測(cè)量中，低通濾波器截止頻率與分辨率的關(guān)系有[22]：

4 實(shí)測(cè)數(shù)據(jù)處理試驗(yàn)及分析

利用航空重力儀在我國(guó)南部某省進(jìn)行了近海重力測(cè)量，本試驗(yàn)所用數(shù)據(jù)為第五架次2號(hào)測(cè)線的測(cè)量數(shù)據(jù)，該測(cè)線是一條東西向的測(cè)線。測(cè)量選取在氣候條件較好的夜間進(jìn)行，以減小氣流對(duì)測(cè)量結(jié)果的影響。測(cè)量過程中，設(shè)計(jì)飛行高度為400 m，飛行速度為60 m/s，實(shí)際飛行測(cè)量高度在設(shè)計(jì)高度的±5米范圍內(nèi)波動(dòng)，圖1為該條測(cè)線的飛行高度變化示意圖。

分別利用卡爾曼平滑濾波器和FIR低通濾波器對(duì)測(cè)量數(shù)據(jù)進(jìn)行了處理，其中在FIR低通濾波器設(shè)計(jì)中采用了海明窗，窗口階數(shù)選取為600，截止頻率選取在0.004 Hz。此外，由于暫時(shí)無法獲得地面重力測(cè)量數(shù)據(jù)，試驗(yàn)過程中以商用軟件的處理結(jié)果作為參考標(biāo)準(zhǔn)，將本試驗(yàn)的處理結(jié)果與之進(jìn)行對(duì)比研究。

根據(jù)本文第2部分和第3部分的方法對(duì)該測(cè)線處理后，所獲重力異常結(jié)果如圖2所示。圖中，橫坐標(biāo)為GPS測(cè)量時(shí)間；縱坐標(biāo)為處理所得重力異常值；實(shí)線為商用軟件處理所得結(jié)果，在本試驗(yàn)中作為參考標(biāo)準(zhǔn)；虛線是FIR低通濾波的重力異常曲線，為方便對(duì)比，受邊界效應(yīng)影響的測(cè)量點(diǎn)已經(jīng)去除；點(diǎn)劃線為卡爾曼平滑濾波處理所得重力異常曲線。

圖1 飛行高度變化示意圖

圖2 FIR低通濾波與卡爾曼平滑濾波處理重力異常對(duì)比圖

表1 FIR低通濾波與卡爾曼平滑濾波處理重力異常誤差對(duì)比表

由圖2可見，經(jīng)過處理后，測(cè)量結(jié)果中的噪聲已經(jīng)被抑制，重力信號(hào)被保留了下來，獲得了高信噪比的重力異常信號(hào)。對(duì)比兩種處理方法的結(jié)果，F(xiàn)IR低通濾波的結(jié)果與參考標(biāo)準(zhǔn)的偏差較為明顯。與此同時(shí)，還應(yīng)注意到無論是FIR處理的結(jié)果還是卡爾曼濾波平滑處理的結(jié)果，重力異常曲線在測(cè)線末端均出現(xiàn)了較為強(qiáng)烈的波動(dòng)，且與參考標(biāo)準(zhǔn)出現(xiàn)了明顯的偏差，分析其原因可能是因?yàn)轱w行過程中產(chǎn)生了較為明顯的垂直方向的機(jī)動(dòng)而引入了部分的噪聲。

表1為兩種濾波方法處理所得重力異常與參考標(biāo)準(zhǔn)的誤差統(tǒng)計(jì)。由表可見，利用FIR低通濾波器與參考標(biāo)準(zhǔn)的均方差為1.2329 mGal，而卡爾曼平滑濾波的處理結(jié)果的均方差僅為0.7219 mGal，遠(yuǎn)遠(yuǎn)優(yōu)于FIR低通濾波的處理效果，這表明相比較于傳統(tǒng)的低通濾波處理方法，卡爾曼平滑濾波處理所獲得的重力異常具有更高的精度。

5 結(jié)論

本文依據(jù)航空重力測(cè)量基本原理，研究了兩種用于垂直加速度改正處理的濾波方法，并對(duì)它們的處理結(jié)果進(jìn)行了對(duì)比試驗(yàn)。由試驗(yàn)結(jié)果得出如下結(jié)論：

（1）本文所提出的卡爾曼平滑濾波狀態(tài)空間模型，經(jīng)實(shí)測(cè)數(shù)據(jù)驗(yàn)證，可以獲取高精度的重力異常信號(hào)，誤差小于1 mGal；

（2）盡管基于低通濾波的處理方法所獲得的重力異常誤差高于卡爾曼平滑濾波處理的結(jié)果，但是其處理方法和步驟相對(duì)簡(jiǎn)單，也不失為一種在實(shí)際應(yīng)用中可以選擇的處理策略；

（3）無論采用何種濾波處理方法，測(cè)量過程中，飛行載體都應(yīng)盡可能保持平穩(wěn)的飛行狀態(tài)，這有助于提高航空重力測(cè)量的精度。

[1] 夏哲仁,石磬,孫中苗.航空重力測(cè)量系統(tǒng)CHAGS[J].測(cè)繪學(xué)報(bào),2004,33(3):216-220.

[2] 夏哲仁,孫中苗.航空重力測(cè)量及其技術(shù)應(yīng)用[J].測(cè)繪科學(xué),2006,31(6):43-46.

[3] 熊盛青.我國(guó)航空重磁勘探技術(shù)現(xiàn)狀與發(fā)展趨勢(shì)[J].地球物理學(xué)進(jìn)展,2009,24(1):113-117.

[4] 熊盛青,周錫華,郭志宏，等.航空重力勘探理論方法及應(yīng)用[M].北京:地質(zhì)出版社,2010.

[5] 黃謨濤,寧津生,歐陽(yáng)永忠，等.航空重力測(cè)量厄特弗斯改正公式注記[J].測(cè)繪學(xué)報(bào),2015,44(1):6-12.

[6] 黃謨濤,寧津生,歐陽(yáng)永忠,等.?？罩亓y(cè)量平臺(tái)傾斜改正模型等價(jià)性證明與驗(yàn)證[J].武漢大學(xué)學(xué)報(bào)（信息科學(xué)版）,2016,41(6):738-744.

[7] 孫中苗, 石磬, 夏哲仁，等. 利用GPS和數(shù)字濾波技術(shù)確定航空重力中的垂直加速度[J]. 測(cè)繪學(xué)報(bào), 2004, 33(2): 110-115.

[8] HAMMADA Y. A comparison of filtering techniques for airborne gravimetry[D]. Calgary: University of Calgary, 1996.

[9] 徐澤揚(yáng). 航空重力數(shù)據(jù)濾波處理方法研究[D].南京:東南大學(xué),2013.

[10] 蔡劭琨. 航空重力矢量測(cè)量及誤差分離方法研究[D].長(zhǎng)沙:國(guó)防科技大學(xué),2014.

[11] 郭志宏, 羅鋒, 王明，等. 航空重力數(shù)據(jù)無限脈沖響應(yīng)低通數(shù)字濾波器設(shè)計(jì)與實(shí)驗(yàn)研究[J]. 地球物理學(xué)報(bào), 2011, 54(8): 2148-2153.

[12] 李宜龍, 曾敏, 孫科，等. 航空重力測(cè)量中FIR低通濾波與高斯濾波的比較[J]. 海洋測(cè)繪, 2014, 34(6):40-42.

[13] 郎駿健，梁星輝，柳林濤，等．航空重力傅里葉基追蹤低通濾波方法研究[J]．地球物理學(xué)報(bào), 2018, 61(12):4737-4745.

[14] BOLOTIN Y V, POPELENSKY M Y. The analysis of accuracy of airborne gravimetry with stochastic models[J]. Aerokosmichekoe priborostroenie, 2003, 4: 42-48.

[15] BOLOTIN Y , DOROSHIN D . Adaptive filtering in airborne gravimetry with hidden Markov chains.[J]. IFAC proceedings volumes,2011,44(1):9996-10001.

[16] BOLOTIN Y V. GOLVAN A A. Methods of inertial gravimetry[J]. Moscow university mechanics bulletin,2013,68(5):117-125.

[17] 鄭崴，張貴賓.自適應(yīng)卡爾曼濾波在航空重力異常解算的應(yīng)用研究[J].地球物理學(xué)報(bào),2016,59(4):275-1283.

[18] 鄭崴,張貴賓,陳濤,等.Sage自適應(yīng)濾波在航空重力數(shù)據(jù)處理的應(yīng)用[J].物探與化探,2015,39(S1):84 -90．

[19] BOLOTIN Y V. Mathematics behind GTGRAV[R]. Moscow: Moscow Lomonosov State University,2009.

[20] 蔡體菁,周薇,鞠玲玲.平臺(tái)式重力儀測(cè)量數(shù)據(jù)的卡爾曼濾波處理[J].中國(guó)慣性技術(shù)學(xué)報(bào),2015,23(6):718-723.

[21] 秦永元,張洪鉞,汪叔華.卡爾曼濾波與組合導(dǎo)航原理（第二版）[M]．西安:西北工業(yè)大學(xué)出版社,2012．

[22] 孫中苗.航空重力測(cè)量理論、方法及應(yīng)用研究[D].鄭州:解放軍信息工程大學(xué), 2004.

[23] 何海波,楊元喜,孫中苗.幾種GPS測(cè)速方法的比較[J].測(cè)繪學(xué)報(bào),2002,31(3):217-221.

[24] 郭志宏,羅鋒,安戰(zhàn)鋒.航空重力數(shù)據(jù)窗函數(shù)法FIR低通數(shù)字濾波試驗(yàn)[J]．物探與化探,2007,31(6):568-571.

Comparing Research for Vertical Acceleration Correction Filtering Methods in Airborne Gravimetry

ZHENG Wei, WANG Jian-ling, CHEN Xue-feng

(School of Electronic Information Engineering, Henan Institute of Technology, Xinxiang 453003, China)

Based on the principle of airborne gravimetry, the vertical acceleration correction methods are demonstrated in this paper. Then, the methods based on low-pass filter and Kalman smoother are designed for the correction respectively. According to the diversities of vertical acceleration correction methods, the filter and the smoother are applied in measurement data processing which is from practical airborne gravimetry. The result of the processing shows that both of the filter and smoother can process the measured data to obtain the gravity anomaly and the errors is about 1 mGal. Specifically, the error based on Kalman smoother is less than 1 mGal.

airborne gravimetry; vertical acceleration correction; low-pass filtering; Kalman smoothing

V241

A

2096–7772(2020)01–0028–05

2019-12-14

鄭崴（1984―），男，河南正陽(yáng)人，講師，博士，主要從事航空重力測(cè)量技術(shù)方法、復(fù)雜信號(hào)處理與提取技術(shù)研究。

（責(zé)任編輯呂春紅）