基于正反Kalman濾波的動(dòng)基座重力儀數(shù)據(jù)處理方法

熊志明，郭 剛，李東明，羅 騁，李海兵

（北京航天控制儀器研究所，北京100039)

基于正反Kalman濾波的動(dòng)基座重力儀數(shù)據(jù)處理方法

熊志明，郭 剛，李東明，羅 騁，李海兵

（北京航天控制儀器研究所，北京100039)

重力數(shù)據(jù)處理是動(dòng)基座重力儀的核心技術(shù)，采用了一種基于正反Kalman濾波的數(shù)據(jù)處理方法提取重力異常值。以動(dòng)基座重力儀（Sea and Air Gravity，SAG)為研究對象，根據(jù)系統(tǒng)參數(shù)推導(dǎo)了Kalman濾波方程，并運(yùn)用正反Kalman濾波方法處理了SAG某飛行架次的數(shù)據(jù)。將提取的重力異常值與同機(jī)搭載飛行的俄羅斯高精度重力儀GT?1A的結(jié)果進(jìn)行比對，試驗(yàn)結(jié)果表明，兩者濾波結(jié)果差值的均方根誤差要小于1mGal。

動(dòng)基座重力儀；重力異常；重力數(shù)據(jù)處理；正反Kalman濾波

0 引言

航空重力測量是在機(jī)載重力和GPS組合系統(tǒng)下進(jìn)行空中重力測量的技術(shù)，它是陸地重力和海洋重力測量的拓展和補(bǔ)充，也是重力測量的重要發(fā)展方向之一。基于測量速度快、覆蓋范圍大等優(yōu)點(diǎn)，航空重力測量可用于大地水準(zhǔn)面測定、無人區(qū)重力測量、遠(yuǎn)程武器發(fā)射區(qū)快速重力測量等方面。目前，航空重力測量技術(shù)在地球物理勘探方面得到廣泛應(yīng)用，服務(wù)領(lǐng)域涉及基礎(chǔ)地質(zhì)調(diào)查與研究、空間科學(xué)、石油、天然氣及固體礦產(chǎn)資源勘探、大地測量等方面［1］。

重力數(shù)據(jù)處理對獲取高精度重力異常值起到了重要作用。航空重力是在飛機(jī)上進(jìn)行測量，飛機(jī)的高頻振動(dòng)將不可避免對重力儀測量數(shù)據(jù)和DGPS解算數(shù)據(jù)產(chǎn)生高頻干擾影響，因此提取的重力異常初值也含有大量的高頻噪聲。為了得到有用的重力異常信號，需進(jìn)行低通濾波消除高頻干擾。本文圍繞重力數(shù)據(jù)低通濾波處理這一核心技術(shù)展開工作。

目前數(shù)據(jù)低通濾波處理方法包括無限脈沖響應(yīng)（Infinite Impulse Response，IIR)低通濾波方法、有限脈沖響應(yīng)（Finite Impulse Response，F(xiàn)IR)低通濾波方法、RC濾波方法、高斯濾波方法、FFT濾波方法等。以上濾波處理方法都是基于非模型的低通濾波算法，容易導(dǎo)致濾波信號中重力異常信號的損失。因此，對于分辨率要求較高的地球物理勘探來說，需選擇基于模型的Kalman濾波估計(jì)方法技術(shù)作進(jìn)一步處理［11］，即對重力異常建模并通過Kalman濾波估計(jì)重力異常值。

本文在文獻(xiàn)［9］和文獻(xiàn)［13］的基礎(chǔ)上對重力異常模型進(jìn)行改進(jìn)，用正反Kalman濾波方法處理SAG的數(shù)據(jù)，并與用標(biāo)準(zhǔn)離散Kalman濾波方法處理的結(jié)果進(jìn)行對比。

1 SAG系統(tǒng)組成及數(shù)據(jù)處理流程

1．1 系統(tǒng)組成

SAG是基于捷聯(lián)慣導(dǎo)的動(dòng)基座重力儀，以捷聯(lián)慣性導(dǎo)航為基礎(chǔ)，結(jié)合差分GPS系統(tǒng)，可以實(shí)現(xiàn)動(dòng)態(tài)重力測量。動(dòng)基座重力測量系統(tǒng)功能組成如圖1所示，分為載體、地面和離線處理工作站3個(gè)組成部分。

圖1 動(dòng)基座重力測量系統(tǒng)功能組成示意圖Fig．1 Schematic diagram of moving base gravity system

載體部分，以動(dòng)基座重力儀為核心單元，同時(shí)還包括由GPS天線和GPS接收機(jī)構(gòu)成的差分GPS移動(dòng)站、顯控記錄單元、電源管理單元。

地面部分，主要是由GPS天線和GPS接收機(jī)構(gòu)成的差分GPS基站。

離線處理工作站部分，主要是對數(shù)據(jù)進(jìn)行離線處理的計(jì)算機(jī)及程序軟件。

SAG系統(tǒng)包括重力儀及后處理軟件系統(tǒng)。重力儀如圖2所示，主要包含慣性測量單元、GPS基站、數(shù)據(jù)采樣記錄儀、28V直流電源、不間斷電源系統(tǒng)UPS、電源顯示控制單元。其中，慣性測量單元如圖3所示。

圖2 重力儀Fig．2 Gravimeter

圖3 慣性測量單元Fig．3 Inertial measurement unit

1．2 數(shù)據(jù)處理流程

數(shù)據(jù)處理流程圖如圖4所示，重力儀輸出的數(shù)據(jù)通過記錄儀存儲，數(shù)據(jù)包括陀螺、加速度計(jì)的原始脈沖以及GPS原始數(shù)據(jù)。GPS原始數(shù)據(jù)經(jīng)差分處理得到精確的位置和速度信息；陀螺和加速度計(jì)的原始脈沖經(jīng)誤差補(bǔ)償后得到加速度和角速度，經(jīng)捷聯(lián)慣性導(dǎo)航得到位置、速度和姿態(tài)。將慣性導(dǎo)航的結(jié)果與GPS的速度和位置進(jìn)行組合導(dǎo)航，來修正姿態(tài)信息得到導(dǎo)航系下的3個(gè)比力值。對重力異常進(jìn)行建模后聯(lián)立導(dǎo)航系下的比力值建立Kalman濾波方程，經(jīng)正反Kalman濾波后得到精確的重力異常值。

由圖4可以看出，數(shù)據(jù)處理流程先進(jìn)行導(dǎo)航解算，再進(jìn)行重力異常值的提取。導(dǎo)航解算處理包括誤差補(bǔ)償、捷聯(lián)慣性解算、組合導(dǎo)航。重力異常提取是在導(dǎo)航解算的基礎(chǔ)下進(jìn)行，包括重力異常建模和Kalman濾波方程構(gòu)建。

圖4 重力數(shù)據(jù)處理流程Fig．4 Flow chart of gravity data processing

2 正反Kalman濾波與標(biāo)準(zhǔn)Kalman濾波比較

2．1 標(biāo)準(zhǔn)離散Kalman濾波

設(shè)離散線性系統(tǒng)的狀態(tài)方程和觀測方程分別為［12］：

其中，Xk為tk時(shí)刻的狀態(tài)向量；?k，k-1為tk-1時(shí)刻至tk時(shí)刻的一步轉(zhuǎn)移陣；Γk-1為系統(tǒng)噪聲耦合矩陣；Zk為tk時(shí)刻的觀測向量；Hk為量測陣；Vk為量測噪聲序列；Wk為系統(tǒng)激勵(lì)噪聲序列。

Kalman濾波要求｛Wk｝和｛Vk｝是互不相關(guān)的零均值的白噪聲序列，有：

Qk和Rk分別為系統(tǒng)噪聲和量測噪聲的方差矩陣，在Kalman濾波中要求它們分別是已知值的非負(fù)定陣和正定陣。

離散Kalman濾波方程為：

其中，Pk為估計(jì)誤差方差陣，Kk為最優(yōu)增益矩陣或者加權(quán)矩陣，Pk／k-1為一步預(yù)測均方誤差，Xk／k-1為狀態(tài)一步預(yù)測。

2．2 正反Kalman濾波器

正反Kalman濾波即將Kalman正向?yàn)V波以及Kalman反向?yàn)V波的結(jié)果進(jìn)行某種形式的加權(quán)求和。將Kalman反向?yàn)V波過程k時(shí)刻所得的一步狀態(tài)預(yù)測量、狀態(tài)估計(jì)分別記為xf（k／k-1)、xf（k)，相應(yīng)的誤差協(xié)方差矩陣分別記為Pf（k／k-1)、Pf（k)。假設(shè)Kalman正向?yàn)V波產(chǎn)生的狀態(tài)估計(jì)量的數(shù)據(jù)個(gè)數(shù)為n，那么Kalman反向?yàn)V波過程的一步狀態(tài)預(yù)測量及其誤差協(xié)方差矩陣初值記為xf（n／n-1)＝0，Pf（n／n-1)＝0。當(dāng)k由n時(shí)刻反向變化為0時(shí)刻時(shí)，正反Kalman反向?yàn)V波的校正過程以及預(yù)測過程表達(dá)式如下所述。

正反Kalman反向?yàn)V波校正過程［9］如式（4)所示。

預(yù)測過程［9］如式（5)所示。

正反Kalman加權(quán)［9］如式（6)所示。

正反Kalman濾波的狀態(tài)量是Kalman正向?yàn)V波的狀態(tài)估計(jì)量xk與Kalman反向?yàn)V波的狀態(tài)估計(jì)量xf（k／k-1)的加權(quán)求和，兩者的權(quán)重分別為

2．3 正反Kalman濾波器的特點(diǎn)

標(biāo)準(zhǔn)離散Kalman濾波是正向Kalman濾波過程，一步狀態(tài)估計(jì)量和觀測量通過Kalman增益以及相應(yīng)的誤差協(xié)方差矩陣各自對狀態(tài)估計(jì)量進(jìn)行修正；正反Kalman濾波器在標(biāo)準(zhǔn)離散Kalman濾波器的基礎(chǔ)上進(jìn)行了一次反向?yàn)V波，并將兩種結(jié)果進(jìn)行加權(quán)求和，這樣能充分使用標(biāo)準(zhǔn)離散Kalman濾波所產(chǎn)生的一步狀態(tài)估計(jì)量、狀態(tài)估計(jì)量及相應(yīng)的誤差協(xié)方差矩陣的信息，可以得到更優(yōu)化的濾波結(jié)果。

3 重力異常建模

SAG是基于捷聯(lián)慣導(dǎo)的動(dòng)基座重力儀，原始輸出量為陀螺儀在載體系下的角度增量和加速度計(jì)在載體系下的速度增量。經(jīng)過SINS／DGPS組合導(dǎo)航濾波方法計(jì)算精確的導(dǎo)航系下的比力f′1、f′2、 f′3。但由于計(jì)算誤差等因素影響，f′3與真值會(huì)有微小偏差，對其建立如下方程：

其中，k1、k2為姿態(tài)失準(zhǔn)角，k3為標(biāo)度因數(shù)誤差，τ3為延時(shí)參數(shù)，δf3為重力噪聲。航空重力標(biāo)量測量的數(shù)學(xué)模型為：

其中，g0為正常重力值，gE是由地球自轉(zhuǎn)和載體相對地球運(yùn)動(dòng)引起的，稱為厄特弗斯改正項(xiàng)。令：

其中，v3為GPS測量的天向速度，Δv3為天向速度誤差。將式（7)和式（9)代入式（8)中，得：

由此建立動(dòng)基座重力儀的高度誤差方程：

其中，Δh為高度差，為GPS測量的高度與IMU積分得到的高度差。Δh′的結(jié)果由式（12)得到：

其中，h′為天向運(yùn)動(dòng)加速度a0經(jīng)過二次積分得到，h″為對GPS天向速度v3積分得到。

重力異常是個(gè)隨機(jī)過程，目前，通常采用高斯?馬爾科夫隨機(jī)過程對重力擾動(dòng)進(jìn)行建模，3階高斯?馬爾科夫過程的微分方程為［10］：

其中，β為相關(guān)參數(shù)，w為輸入的白噪聲。

根據(jù)式（13)，重力異常的模型可以寫成如下形式：

qg是方差強(qiáng)度為Qg的白噪聲，Ag、Bg、Cg為常值矩陣。

式（14)結(jié)合動(dòng)基座重力儀的高度誤差方程，即可得到重力異常值估計(jì)的Kalman濾波方程：

式中，狀態(tài)變量為x＝［xaxg］T，xa＝［Δv3Δh τ3k1k2k3］T，qf、qk、δh為白噪聲。高度差Δh′為狀態(tài)觀測量。

4 飛行試驗(yàn)結(jié)果對比

采用上述重力異常模型，分別用標(biāo)準(zhǔn)Kalman濾波和正反Kalman濾波方法對SAG某次飛行試驗(yàn)數(shù)據(jù)進(jìn)行了處理，本文給出了其中6條飛行測線處理結(jié)果，并將兩種結(jié)果與同機(jī)飛行的GT重力儀結(jié)果進(jìn)行對比。飛機(jī)平穩(wěn)等高度飛行，飛行的平均高度為1500m，飛行平均速度為65m／s。飛行軌跡如圖5所示。

圖5 飛行軌跡圖Fig．5 Flight path

由圖5可以看出，飛行軌跡為6條南北測線。從左至右6條測線分別為測線1至測線6。6條測線的對比曲線分別如圖6～圖11所示。

圖6 測線1對比曲線Fig．6 Comparison of line 1

圖7 測線2對比曲線Fig．7 Comparison of line 2

圖8 測線3對比曲線Fig．8 Comparison of line 3

圖6～圖11中的虛線為采用本文的重力異常模型和正反Kalman濾波方法處理的SAG的重力異常曲線（SAGFB)；點(diǎn)線為用標(biāo)準(zhǔn)Kalman濾波方法處理SAG的重力異曲線（SAGK)；實(shí)線為GT?1A的重力異常曲線（GT?1A)；曲線的橫坐標(biāo)為UTC時(shí)間，單位為s；縱坐標(biāo)為重力異常值，單位為mGal（10-5m／s2)。由曲線可以看出，兩種Kalman濾波方法估計(jì)出來的SAG的重力異常曲線與GT?1A的重力異常曲線趨勢基本一致。為了更好地比較兩種濾波方法，本文統(tǒng)計(jì)了兩種濾波結(jié)果與GT?1A結(jié)果的符合性精度，如表1所示。

圖9 測線4對比曲線Fig．9 Comparison of line 4

圖10 測線5對比曲線Fig．10 Comparison of line 5

表1 SAG與GT?1A測量結(jié)果的符合性統(tǒng)計(jì)表Table 1 Statistical table of measuring result between SAG and GT?1A

表1以SAG與GT?1A在每條測線上所有測點(diǎn)的均方根誤差來評價(jià)2套不同重力儀的符合性，并比較兩種濾波方法的精度，測點(diǎn)的采樣周期為2Hz。每條測線的測點(diǎn)數(shù)如表1中第6列所示，各條測線差值的均方根在表的第3列和第5列給出。由表1可以看出，用正反Kalman濾波方法處理的各條測線的差值均方根在1mGal以內(nèi)，6條測線誤差的均方根的平均值為0.9073mGal；用標(biāo)準(zhǔn)離散Kalman濾波方法處理的各條測線的差值均方根高于1mGal，均方根的平均值為1.1273mGal。由此看出，正反Kalman濾波的結(jié)果要好于標(biāo)準(zhǔn)離散Kalman濾波的結(jié)果。

每條測線的誤差均值A(chǔ)計(jì)算公式如式（16)所示。

其中，n為測線的點(diǎn)個(gè)數(shù)，即表格第4列所示；x1、x2、…、xn為GT?1A處理的各個(gè)測點(diǎn)的重力異常值；y1、y2、…、yn為SAG處理的各個(gè)測點(diǎn)的重力異常值；a1、a2、…、an為兩種儀器處理的各個(gè)測點(diǎn)的重力異常之差。

誤差的均方根D的計(jì)算公式如式（17)所示。

5 結(jié)論

本文探討了基于正反Kalman濾波的重力異常提取方法，采用該方法對動(dòng)基座重力儀SAG實(shí)際飛行測量數(shù)據(jù)進(jìn)行了處理，并與同機(jī)飛行的國外進(jìn)口GT?1A重力儀結(jié)果進(jìn)行了對比。根據(jù)SAG與GT?1A的符合性統(tǒng)計(jì)結(jié)果表明，本文的濾波方法獲得的處理精度與GT?1A精度相當(dāng)，能滿足航空重力勘察測繪的要求。這驗(yàn)證了重力異常模型的合理性和正反Kalman濾波算法的實(shí)用性，也從側(cè)面反映了SAG已達(dá)到國外進(jìn)口產(chǎn)品同等精度。下一步研究的重點(diǎn)在于進(jìn)一步優(yōu)化重力異常模型，得到更精確的模型。

本次飛行試驗(yàn)得到了中國國土資源部物探遙感中心的大力幫助，在此表示誠摯的感謝。

［1］熊盛青，周錫華，郭志宏，等.航空重力勘探理論方法及應(yīng)用［M］.北京：地質(zhì)出版社，2010. XIONG Sheng?qing，ZHOU Xi?hua，GUO Zhi?hong，et al.Thoery，method and application of the airborne gravity prospecting［M］.Beijing：Geology Press，2010.

［2］郭志宏，羅鋒，王明，等.航空重力數(shù)據(jù)無限脈沖相應(yīng)低通數(shù)字濾波器設(shè)計(jì)與實(shí)驗(yàn)研究［J］.地球物理學(xué)報(bào)，2011，54（8）：2148?2153. GUO Zhi?hong，LUO Feng，WANG Ming，et al.The de?sign and experiment of IIR lowpass digital filters for air?borne gravity data［J］.Chinese Journal of Geophysics，2011，54（8）：2148?2153.

［3］蔡體菁，周薇，丁昊，等.航空重力數(shù)據(jù)與重力梯度信號處理［C］.慣性技術(shù)發(fā)展動(dòng)態(tài)發(fā)展方向研討會(huì)，2014：53?56. CAI Ti?jing，ZHOU Wei，DING Hao，et al.Data pro?cessing for airborne gravity and land gravity gradient［C］. Dynamic Development Direction Seminar of Inertial Tech?nology development，2014：53?56.

［4］田顏鋒，李?yuàn)檴?，劉曉剛，?航空重力數(shù)據(jù)低通濾波處理及其功效分析［J］.大地測量與地球動(dòng)力學(xué)，2009，29（6）：130?132. TIAN Yan?feng，LI Shan?shan，LIU Xiao?gang，et al. Low pass filter and its effect in data processing of aerial gravity measurements［J］.Journal of Geodesy and Geody?namics，2009，29（6）：130?132.

［5］郭志宏，羅鋒，安戰(zhàn)鋒.航空重力數(shù)據(jù)窗函數(shù)法FIR低通數(shù)字濾波試驗(yàn)［J］.物探與化探，2007，31（6）：568?571. GUO Zhi?hong，LUO Feng，AN Zhan?feng.Experiment researches on FIR low pass digital filters based on window functions of airborne gravity data［J］.Geophysical And Geochemical Exploration，2007，31（6）：568?571.

［6］李宜龍，曾敏，孫科，等.航空重力測量中FIR濾波與高斯濾波的比較［J］.海洋測繪，2014，34（6）：40?42. LI Yi?long，ZENG Min，SUN Ke，et al.Comparison be?tween FIR filter and Gaussian filter in airborne gravity sur?veys［J］.Hydrographic Surveying and Charting，2014，34（6）：40?42.

［7］袁園.航空重力數(shù)據(jù)常用濾波方法對比研究［J］.吉林大學(xué)學(xué)報(bào)，2010，40（S1）：6?9＋38. YUAN Yuan.Comparative study of common filtering methods used in airborne gravity data［J］.Journal of Jilin University，2010，40（S1）：6?9＋38.

［8］孫中苗，夏哲仁.FIR低通差分器的設(shè)計(jì)及其在航空重力測量中的應(yīng)用［J］.地球物理學(xué)報(bào)，2000，43（6）：850?855. SUN Zhong?miao，XIA Zhe?ren.Design of FIR lowpass differentiator and its applications in airborne gravimetry［J］. Chinese Journal of Geophysics，2000，43（6）：850?855.

［9］孟濤.卡爾曼濾波在航空重力數(shù)據(jù)處理中的應(yīng)用研究［D］.東南大學(xué)，2012. MENG Tao.Kalman filter for airborne gravity data［D］. Southeast University，2012.

［10］Gelb A.Applied optimal estimation［M］.Cambridge，Massachusetts：MIT Press，1974.

［11］王靜波，熊盛青，郭志宏，等.航空重力數(shù)據(jù)Kalman濾波平滑技術(shù)應(yīng)用研究［J］.地球物理學(xué)進(jìn)展，2012，27（4）：1717?1722. WANG Jing?bo，XIONG Sheng?qing，GUO Zhi?hong，et al.Kalman smoothing for airborne gravity data［J］.Pro?gress in Geophysics，2012，27（4）：1717?1722.

［12］秦永元，張洪鉞，汪叔華.卡爾曼濾波與組合導(dǎo)航原理［M］.西安：西北工業(yè)大學(xué)出版社，1998. QIN Yong?yuan，ZHANG Hong?yue，WANG Shu?hua.Kal?man filter and principal of integrated navigation［M］.Xi'an：North Western Polytechnical University Press，1998.

［13］王靜波，熊盛青，郭志宏，等.航空重力異常估計(jì)方法研究［J］.物探與化探，2011，35（4）：493?498. WANG Jing?bo，XIONG Sheng?qing，GUO Zhi?hong，et al.Research on methods for estimating airborne gravity anomalies［J］.Geophysical＆Geochemical Exploration，2011，35（4）：493?498.

［14］Senobari M S.New results in airborne vector gravimetry using strapdown INS／DGPS［J］.Journal of Geodesy，2010，84（5）：277?291.

［15］李東明，郭剛，薛正兵，等.激光捷聯(lián)慣導(dǎo)車載重力測量試驗(yàn)［J］.導(dǎo)航與控制，2013，12（4）：75?78. LI Dong?ming，GUO Gang，XUE Zheng?bing，et al. Gravimetry test of laser strapdown INS on ground?vehicle［J］.Navigation and Control，2013，12（4）：75?78.

［16］孫中苗，翟振和，李迎春.航空重力儀發(fā)展現(xiàn)狀和趨勢［J］.地球物理學(xué)進(jìn)展，2013，28（2）：1?8. SUN Zhong?miao，ZHAI Zhen?he，LI Ying?chun.Status and development of airborne gravimeter［J］.Progress in Geophysics，2013，28（2）：1?8.

Gravity Data Processing Based on Forward and Backward Kalman Filter

XIONG Zhi?ming，GUO Gang，LI Dong?ming，LUO Cheng，LI Hai?bing
（Beijing Institute of Aerospace Control Devices，Beijing 100039)

Gravity data processing is the key technology of moving base gravimeter.A data processing method based on forward and backward Kalman filter was used to extract gravity anomaly.According to the system parameter of sea and air gravity（SAG)，Kalman filter equation was deduced and flight data of SAG was processed through forward and backward Kalman filter.It was shown that the root?mean?square error of the difference of gravity anomaly is lower than 1 mGal by comparing the gravity anomaly result of SAG with GT?1A which both been carried on the airplane.

moving base gravimeter;gravity anomaly;gravity data processing;forward and backward Kalman filter

P631.1＋25

A

1674?5558（2017)03?01232

10.3969／j.issn.1674?5558.2017.01.014

熊志明，男，碩士，研究方向?yàn)閯?dòng)基座重力測量系統(tǒng)。

2016?01?13

國家高技術(shù)研究發(fā)展（863)計(jì)劃（編號：2011AA060506)；國家國際科技合作專項(xiàng)（編號：2014DFR80750)；航天科技集團(tuán)公司九院創(chuàng)新基金項(xiàng)目（動(dòng)基座重力測量系統(tǒng)，航空重磁一體化綜合信息系統(tǒng))