鄭崴， 張貴賓

1 中國(guó)地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院， 北京　100083 2 河南科技學(xué)院信息工程學(xué)院， 河南新鄉(xiāng)　453003

自適應(yīng)卡爾曼濾波在航空重力異常解算的應(yīng)用研究

鄭崴1,2， 張貴賓1*

1 中國(guó)地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院， 北京100083 2 河南科技學(xué)院信息工程學(xué)院， 河南新鄉(xiāng)453003

摘要依據(jù)航空重力測(cè)量基本原理，構(gòu)建了航空重力異常解算的卡爾曼濾波模型，將新息自適應(yīng)卡爾曼濾波器(IAE，Innovation based Adaptive Estimation)應(yīng)用于量測(cè)噪聲未知的航空重力異常解算.針對(duì)IAE濾波器滑動(dòng)窗口寬度難以準(zhǔn)確確定的問(wèn)題，通過(guò)對(duì)多個(gè)不同滑動(dòng)窗口新息協(xié)方差估計(jì)的加權(quán)平均，獲得改進(jìn)的IAE濾波器，該IAE濾波器不僅具有量測(cè)噪聲自適應(yīng)估計(jì)能力，還能實(shí)現(xiàn)滑動(dòng)采樣窗口的優(yōu)化選取.試驗(yàn)結(jié)果表明，IAE濾波器可以降低因量測(cè)噪聲統(tǒng)計(jì)信息不明引起的解算誤差，改進(jìn)IAE解算的重力異常誤差約為1 mGal.

關(guān)鍵詞航空重力異常解算； 自適應(yīng)卡爾曼濾波； 新息； 量測(cè)噪聲

1引言

航空重力測(cè)量是以飛機(jī)為載體，測(cè)量近地空中重力異常的重力測(cè)量技術(shù)，它綜合運(yùn)用了航空重力儀，全球定位系統(tǒng)GPS(Global Position System)，慣性導(dǎo)航系統(tǒng)INS(Inertial Navigation System)，無(wú)線電等技術(shù)設(shè)備(張昌達(dá)，2005；熊盛青，2009a；王靜波等，2009，2010).航空重力測(cè)量數(shù)據(jù)經(jīng)過(guò)后處理，便可得到所需的航空重力異常，這個(gè)數(shù)據(jù)后處理的過(guò)程稱為航空重力異常解算(簡(jiǎn)稱為航重解算).近年來(lái)，差分GPS技術(shù)的進(jìn)步使得航空重力測(cè)量精度大幅度提高，解算所得重力異常精度大部分可達(dá)到1～2 mGal(1 mGal=10-5m·s-2)，最高的已達(dá)0.5 mGal(郭志宏等，2008；熊盛青，2009b; 熊盛青等，2010；Cai et al.，2013).

航空重力測(cè)量是在高速飛行運(yùn)動(dòng)中進(jìn)行的，測(cè)量數(shù)據(jù)必然會(huì)受到包括飛機(jī)運(yùn)動(dòng)產(chǎn)生的擾動(dòng)加速度的干擾，干擾噪聲的幅度最高可達(dá)到106mGal以上(孫中苗，2004)，遠(yuǎn)遠(yuǎn)高于幾十到幾百毫伽的重力異常信號(hào).航空重力測(cè)量數(shù)據(jù)的噪聲主要體現(xiàn)為高頻干擾，依據(jù)測(cè)量數(shù)據(jù)的頻域特點(diǎn)設(shè)計(jì)低通濾波器可消除噪聲影響提取所需的重力異常信號(hào).目前，能夠?qū)崿F(xiàn)線性相位的FIR(Finite Impulse Response)濾波器是航重解算普遍采用的低通濾波器，利用FIR低通濾波器可以獲得精度約1 mGal左右的航空重力異常(郭志宏等，2007, 2011).FIR濾波器的顯著優(yōu)點(diǎn)是具有成熟的設(shè)計(jì)和實(shí)現(xiàn)方法，但由于重力信號(hào)與噪聲在頻譜上沒(méi)有明顯的分界線，為獲得高精度的解算輸出，需設(shè)計(jì)高階的低通濾波器以分離重力異常信號(hào)與噪聲(蔡劭琨等，2010；Cai et al.，2010；羅峰等，2012)，這加重了邊界效應(yīng)對(duì)解算結(jié)果的影響.

部分學(xué)者致力于利用卡爾曼濾波器代替數(shù)字低通濾波器實(shí)現(xiàn)航空重力異常的解算(Bolotin et al.，2003；Bolotin and Popelensky，2007；Bolotin and Golovan，2013).卡爾曼濾波器通常是采用狀態(tài)空間模型在時(shí)間域?qū)ο到y(tǒng)的運(yùn)動(dòng)狀態(tài)進(jìn)行描述，因此卡爾曼濾波器是基于系統(tǒng)運(yùn)動(dòng)模型的濾波方法.狀態(tài)空間模型的應(yīng)用將不可觀測(cè)的變量(狀態(tài)變量)并入可觀測(cè)的模型，并與其一起處理得到相應(yīng)的估計(jì)輸出，所以卡爾曼濾波可實(shí)現(xiàn)多維信號(hào)的并行處理(王靜波，2010；王靜波等，2011).在進(jìn)行航重解算時(shí)，還能夠?qū)娇罩亓x校準(zhǔn)參數(shù)進(jìn)行實(shí)時(shí)估計(jì)(Zhou and Cai，2013；Zheng et al.，2014).此外，航重解算是在線下進(jìn)行的，可對(duì)卡爾曼濾波的估計(jì)結(jié)果進(jìn)行平滑處理，獲得高精度的航空重力異常輸出(潘炎冰，2012；李瑞，2014).但是，卡爾曼濾波能夠獲得準(zhǔn)確解算輸出的重要前提，是能夠掌握系統(tǒng)噪聲和量測(cè)噪聲的統(tǒng)計(jì)信息，而這一前提往往是難以滿足的.解決這一問(wèn)題的辦法就是用自適應(yīng)卡爾曼濾波器(AKF, Adaptive Kalman filter)代替?zhèn)鹘y(tǒng)的標(biāo)準(zhǔn)卡爾曼濾波器(SK, Standard Kalman filter)進(jìn)行航重解算.基于新息協(xié)方差估計(jì)的自適應(yīng)卡爾曼濾波器(IAE, Innovation-based Adaptive Estimation)已證明可以提高噪聲統(tǒng)計(jì)特性不明時(shí)的估計(jì)精度(Mohamed，1999；Mohamed and Schwarz，1999)，IAE的優(yōu)點(diǎn)在于實(shí)現(xiàn)方法簡(jiǎn)單且能夠?qū)α繙y(cè)噪聲協(xié)方差進(jìn)行實(shí)時(shí)估計(jì)，但其新息協(xié)方差估計(jì)器的滑動(dòng)采樣窗口寬度缺乏明確的確定方法.

航重解算中，卡爾曼濾波所需的量測(cè)量由GPS測(cè)定的位置或速度數(shù)據(jù)提供.由于測(cè)量過(guò)程中飛機(jī)運(yùn)動(dòng)狀態(tài)變化劇烈且具有不規(guī)律性，量測(cè)噪聲參數(shù)難以準(zhǔn)確預(yù)設(shè)，因此這里著重討論利用IAE濾波器實(shí)現(xiàn)量測(cè)噪聲自適應(yīng)的航重解算方法.本文將依據(jù)航空重力測(cè)量的基本原理，建立航空重力異常解算的卡爾曼濾波模型，設(shè)計(jì)適用于G型航空重力測(cè)量系統(tǒng)的IAE濾波器.另外，針對(duì)IAE濾波器的滑動(dòng)窗口寬度難以準(zhǔn)確確定的問(wèn)題，采用多個(gè)不同滑動(dòng)窗口寬度的新息協(xié)方差估計(jì)器，通過(guò)對(duì)多個(gè)新息協(xié)方差估計(jì)的加權(quán)平均，獲得滑動(dòng)窗口寬度優(yōu)化的改進(jìn)IAE濾波器.最后，將IAE濾波器和改進(jìn)的IAE濾波器用于實(shí)測(cè)數(shù)據(jù)的航重解算試驗(yàn).

2新息自適應(yīng)卡爾曼濾波器及改進(jìn)

2.1新息自適應(yīng)卡爾曼濾波器

設(shè)在時(shí)刻k，離散線性系統(tǒng)的狀態(tài)方程和量測(cè)方程可分別表示為

(1)

(2)

Xk為k時(shí)刻的狀態(tài)變量；Φk,k-1表示由k-1時(shí)刻到k時(shí)刻的一步轉(zhuǎn)移矩陣；Γk-1為系統(tǒng)噪聲驅(qū)動(dòng)陣；Zk為時(shí)刻k的量測(cè)量；Hk為量測(cè)陣；Wk-1為系統(tǒng)激勵(lì)噪聲序列；Vk為量測(cè)噪聲序列.Qk是系統(tǒng)噪聲的協(xié)方差陣；Rk為量測(cè)噪聲協(xié)方差陣，則標(biāo)準(zhǔn)卡爾曼濾波的算法為(秦永元等，2012)

狀態(tài)一步預(yù)測(cè)：

(3)

一步預(yù)測(cè)誤差：

(4)

濾波增益：

(5)

狀態(tài)估計(jì)：

(6)

估計(jì)均方誤差：

(7)

新息(Innovation)定義為k時(shí)刻濾波器的量測(cè)估計(jì)值與實(shí)際量測(cè)值之差，它表示由量測(cè)量帶來(lái)的、濾波器無(wú)法預(yù)測(cè)的那部分信息，即

(8)

新息序列的理論協(xié)方差為Crk，則有

(9)

當(dāng)系統(tǒng)噪聲和量測(cè)噪聲服從不相關(guān)的高斯分布且新息序列滿足各態(tài)遍歷性時(shí)，可根據(jù)新息序列的滑動(dòng)平均獲得新息協(xié)方差的極大似然估計(jì)，有(MohamedandSchwarz，1999；岳曉奎和袁建平，2005；卞鴻巍等，2006)

(10)

(11)

則量測(cè)噪聲協(xié)方差陣的自適應(yīng)估計(jì)為

(12)

2.2滑動(dòng)采樣窗口寬度的優(yōu)化

實(shí)際應(yīng)用中，如建立的系統(tǒng)模型不準(zhǔn)確，不能真實(shí)反映物理過(guò)程時(shí)，易導(dǎo)致模型與獲得的量測(cè)值不匹配，將會(huì)出現(xiàn)濾波發(fā)散，這時(shí)需通過(guò)增大當(dāng)前時(shí)刻量測(cè)信息在狀態(tài)估計(jì)中的加權(quán)系數(shù)來(lái)抑制濾波發(fā)散(秦永元等，2012).

圖1　滑動(dòng)采樣窗口優(yōu)化Fig.1　Structure of optimization for sliding sampling window

(13)

(14)

(15)

可見(jiàn)，相比于普通的IAE濾波器，改進(jìn)的IAE濾波器配備有多個(gè)不同滑動(dòng)采樣窗口寬度的新息協(xié)方差估計(jì)器，通過(guò)對(duì)不同采樣窗口的新息協(xié)方差估計(jì)進(jìn)行加權(quán)平均，允許新息協(xié)方差估計(jì)器的滑動(dòng)采樣窗口在一定范圍內(nèi)動(dòng)態(tài)變化，一方面可以減弱不當(dāng)?shù)幕瑒?dòng)采樣窗口造成的解算誤差，另一方面也可以提高濾波器工作的穩(wěn)定性，其估計(jì)精度劣于取最優(yōu)滑動(dòng)采樣窗口的情況，但優(yōu)于其他滑動(dòng)采樣窗口寬度的估計(jì)結(jié)果.多個(gè)新息方差估計(jì)器的結(jié)構(gòu)如圖1所示.

2.3 RTS(Rauch-Tung-Striebel)平滑器

(16)

(17)

(18)

3航空重力異常解算的卡爾曼濾波模型

G型航空重力測(cè)量系統(tǒng)是三軸穩(wěn)定平臺(tái)型測(cè)量系統(tǒng)，它包括一個(gè)具有水平穩(wěn)定作用的慣性平臺(tái)(INS)，一個(gè)重力傳感器(航空重力儀)，以及機(jī)載GPS接收器(移動(dòng)站)和地面GPS接收器(接收站).本部分將根據(jù)航空重力儀運(yùn)動(dòng)方程及GPS測(cè)量方程構(gòu)建航重解算的卡爾曼濾波模型.

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

航空重力測(cè)量的基本原理是從航空重力儀測(cè)得的比力數(shù)據(jù)(飛機(jī)運(yùn)動(dòng)加速度與重力加速度之和)中扣除飛機(jī)運(yùn)動(dòng)加速度的影響，獲得所需要的重力加速度(熊盛青等，2010).根據(jù)牛頓第二定律，航空重力儀(垂直加速度計(jì))傳感器探測(cè)質(zhì)量(SM，Sensitive Mass)在地理垂直方向上的動(dòng)力學(xué)方程可表示為(Bolotin，2009)

(19)

厄特缶斯改正是為了補(bǔ)償因地球自轉(zhuǎn)產(chǎn)生的離心力和飛機(jī)運(yùn)動(dòng)產(chǎn)生的離心力附加在重力儀測(cè)量質(zhì)點(diǎn)造成的測(cè)量誤差，可計(jì)算為

(20)

uE為地球自轉(zhuǎn)速度；φ為飛機(jī)所處地理緯度；vE、vN分別是飛機(jī)的東向速度和北向速度；RE、RN分別是地理東向和北向的曲率半徑.

重力正常場(chǎng)和高度改正項(xiàng)G0(φ,h)可由Helmert公式近似計(jì)算(Bolotin，2009)：

(21)

式中，G=9.7803 m·s-2，β=0.005302，

β1=-0.000007，β0=-0.00014 m·s-2，

ω0=0.00123 s-1.

不難看出，為獲得式(19)中的重力異常Δg，需要對(duì)航空重力儀測(cè)量數(shù)據(jù)進(jìn)行一系列的改正.其中，一部分改正項(xiàng)有嚴(yán)密的解析公式，可以直接計(jì)算得到，如厄特缶斯改正、正常場(chǎng)和高度改正，它們統(tǒng)稱為解析類改正；而另一類改正無(wú)法直接用解析公式來(lái)計(jì)算，需借助濾波方法來(lái)實(shí)現(xiàn)，即垂直加速度改正.因此，利用卡爾曼濾波實(shí)現(xiàn)重力異常解算實(shí)際上是在消除噪聲干擾的同時(shí)實(shí)現(xiàn)垂直加速度改正的過(guò)程.3.2卡爾曼濾波模型的建立

將航空重力儀SM的實(shí)時(shí)高程，作為卡爾曼濾波的量測(cè)方程，寫(xiě)為

(22)

h′為重力儀SM所處的實(shí)時(shí)高度；h為SM的實(shí)際高度；δh為測(cè)量誤差.

航空重力儀及兩個(gè)水平加速度計(jì)的測(cè)量方程可表示為(Bolotin，2009)

(23)

將式(23)代入式(19)聯(lián)立并整理得

(24)

其中

(25)

fΣ為已知的作用在SM上的合力產(chǎn)生的加速度總和；qΣ為系統(tǒng)噪聲項(xiàng).為實(shí)現(xiàn)航重解算，將重力異?？醋餮販y(cè)線分布的隨機(jī)過(guò)程，并用濾波器表示，有

(26)

將式(24)和式(26)合并作為卡爾曼濾波的系統(tǒng)方程，式(22)為卡爾曼濾波的量測(cè)方程，并用狀態(tài)空間形式表示為

(27)

(28)

其中，Xf為狀態(tài)向量，重力異常Δg為其分量，即

(29)

Af為動(dòng)力學(xué)矩陣；Bf為系統(tǒng)噪聲驅(qū)動(dòng)矩陣；Cf為控制項(xiàng)驅(qū)動(dòng)矩陣；qf和δh分別為系統(tǒng)噪聲和量測(cè)噪聲項(xiàng)；uf為系統(tǒng)的固定控制項(xiàng).將式(27)和(28)進(jìn)行離散化，再利用設(shè)計(jì)的IAE濾波器和RTS平滑器進(jìn)行處理，即可實(shí)現(xiàn)重力異常的估計(jì).

4航空重力異常解算試驗(yàn)

4.1測(cè)量基本情況

我們?cè)?007年利用G型航空重力測(cè)量系統(tǒng)在某地進(jìn)行了實(shí)地測(cè)量.實(shí)驗(yàn)所用數(shù)據(jù)為第五架次飛行測(cè)量數(shù)據(jù)，飛行時(shí)間為4 h.此次測(cè)量是為了驗(yàn)證系統(tǒng)工作的穩(wěn)定性及重復(fù)測(cè)量結(jié)果的一致性.圖2為測(cè)線位置示意圖，圖中橫軸為相對(duì)經(jīng)度，縱軸為相對(duì)緯度.由圖可知，此次測(cè)量共有6條東西方向的飛行測(cè)線，其中1號(hào)、2號(hào)、3號(hào)及6號(hào)測(cè)線為四條短測(cè)線，它們是沿同一位置往返飛行4次的重復(fù)測(cè)線；4號(hào)測(cè)線和5號(hào)測(cè)線為兩條長(zhǎng)飛行測(cè)線，其中5號(hào)測(cè)線相較其他測(cè)線有約0.002°的北向偏移.測(cè)量選取在氣象條件較好的夜間進(jìn)行，以減小氣流對(duì)測(cè)量結(jié)果的影響.測(cè)量中，設(shè)計(jì)飛行高度為400 m，實(shí)際飛行測(cè)量高度在設(shè)計(jì)飛行高度的±15 m范圍內(nèi)波動(dòng)，圖3為飛行測(cè)量高度變化示意圖.此外，飛行測(cè)量中的設(shè)計(jì)飛行速度為60 m·s-1.

圖2　測(cè)線位置示意圖Fig.2　Position diagram of survey lines

圖3　飛行高度示意圖Fig.3　Chart of flight height

剔除實(shí)測(cè)數(shù)據(jù)中的非測(cè)線數(shù)據(jù)，并進(jìn)行解析類改正后，所得結(jié)果(未濾波)如圖4所示.不難看出，航空重力測(cè)量信號(hào)受噪聲干擾非常嚴(yán)重，幾十到幾百毫伽的重力異常完全被數(shù)萬(wàn)毫伽的噪聲所淹沒(méi).

關(guān)于解算試驗(yàn)有以下兩點(diǎn)需要說(shuō)明，第一，由于暫無(wú)法獲得實(shí)際地面重力測(cè)量數(shù)據(jù)，試驗(yàn)過(guò)程中將航空重力測(cè)量系統(tǒng)配套軟件解算結(jié)果作為解算參考標(biāo)準(zhǔn)，將本試驗(yàn)的解算結(jié)果與參考標(biāo)準(zhǔn)進(jìn)行對(duì)比研究；第二，由于RTS平滑后的重力異常精度高于只使用卡爾曼濾波器估計(jì)的重力異常精度(潘炎冰，2012；李瑞，2014)，因此后文解算結(jié)果如無(wú)特別說(shuō)明均是指平滑處理后的重力異常.

4.2不同預(yù)設(shè)量測(cè)噪聲的航空重力異常解算試驗(yàn)

將SK解算的重力異常與參考標(biāo)準(zhǔn)進(jìn)行對(duì)比，通過(guò)反復(fù)調(diào)整量測(cè)噪聲和系統(tǒng)噪聲，獲得與解算參考標(biāo)準(zhǔn)誤差最小時(shí)的量測(cè)噪聲和系統(tǒng)噪聲分別為R0=(0.0153)2和Q0=diag[0,(3.5×10-5)2,0,0,0,0,0,(0.85×10-7)2]，diag[·]為取對(duì)角線元素符號(hào).實(shí)驗(yàn)中SK和IAE的系統(tǒng)噪聲設(shè)定為Q0，初始量測(cè)噪聲R分別預(yù)設(shè)為R/R0=0.01、0.25、25、100進(jìn)行解算.同時(shí)，為確保新息協(xié)方差估計(jì)的準(zhǔn)確性，IAE的滑動(dòng)采樣窗口寬度暫設(shè)為80.

SK和IAE在不同預(yù)設(shè)量測(cè)噪聲時(shí)，解算的重力異常如圖5所示.為方便分析，將解算所得6條飛行測(cè)線的重力異常按時(shí)間順序連接組成一條曲線，紅色曲線為解算參考標(biāo)準(zhǔn)，藍(lán)色曲線為SK解算的重力異常曲線，黑色曲線為IAE解算的重力異常曲線.由圖可見(jiàn)，SK解算的重力異常曲線隨著預(yù)設(shè)R的變化出現(xiàn)較大的波動(dòng)，而且預(yù)設(shè)的R與R0的偏差越大，SK解算重力異常的誤差也越高；而在同樣的量測(cè)噪聲條件下，IAE解算的重力異常曲線基本保持穩(wěn)定，即IAE可降低因預(yù)設(shè)量測(cè)噪聲不當(dāng)引起的解算誤差.

統(tǒng)計(jì)SK和IAE解算輸出與參考標(biāo)準(zhǔn)的誤差，如表1所示.可以看到，在預(yù)設(shè)量測(cè)噪聲R=0.01R0時(shí)，SK解算結(jié)果的均方誤差高達(dá)3.423367mGal；當(dāng)預(yù)設(shè)量測(cè)噪聲R=0.25R0時(shí)，SK解算重力異常與參考標(biāo)準(zhǔn)的均方誤差為1.198765mGal.在IAE解算誤差一欄中，不同預(yù)設(shè)量測(cè)噪聲條件下，IAE解算重力異常的最大均方誤差為0.923542mGal，最小均方誤差為0.887827mGal，兩者的差距僅為0.035715mGal，而且IAE解算結(jié)果與解算參考標(biāo)準(zhǔn)的誤差都保持在1 mGal以下.所以，IAE濾波器在解算過(guò)程中對(duì)量測(cè)噪聲協(xié)方差進(jìn)行實(shí)時(shí)的估計(jì)和修正，減弱了解算精度對(duì)先驗(yàn)量測(cè)噪聲統(tǒng)計(jì)信息的依賴.因此，相比基于SK的航重解算方法，利用IAE進(jìn)行航重解算更具有普遍性的意義.

圖4　未濾波航空重力異常Fig.4　Un-filtered airborne gravity anomaly

解算濾波方法比較點(diǎn)數(shù)誤差統(tǒng)計(jì)項(xiàng)預(yù)設(shè)量測(cè)噪聲(R/R0)0.010.2525100SK19518平均誤差(mGal)3.0827970.9786831.9484592.57978619518均方誤差(mGal)3.4233671.1987652.0226482.652140IAE(N=80)19518平均誤差(mGal)0.6822390.6834810.7147730.71157719518均方誤差(mGal)0.8878270.8891970.9235420.920205

圖5　SK、IAE與解算參考標(biāo)準(zhǔn)解對(duì)比圖(a) R/R0=0.01； (b) R/R0=0.25； (c) R/R0=25； (d) R/R0=100.Fig.5　Airborne gravity anomaly of SK、IAE VS reference

4.3不同滑動(dòng)窗口的航空重力異常解算試驗(yàn)

由IAE算法步驟可知，滑動(dòng)采樣窗口寬度N直接關(guān)系到新息協(xié)方差估計(jì)的結(jié)果，不當(dāng)?shù)幕瑒?dòng)采樣窗口也會(huì)導(dǎo)致IAE解算精度的下降.本部分將利用不同滑動(dòng)窗口的IAE濾波器和改進(jìn)IAE濾波器進(jìn)行解算試驗(yàn)，驗(yàn)證改進(jìn)IAE濾波器的優(yōu)化效果，并討論N對(duì)IAE解算重力異常精度的影響.試驗(yàn)中，IAE濾波器的滑動(dòng)采樣窗口取N=40、80、120、160、200分別進(jìn)行解算；改進(jìn)IAE濾波器的新息協(xié)方差估計(jì)器設(shè)定為5個(gè)，各新息協(xié)方差估計(jì)器的采樣窗口分別為N1=40，N2=80，N3=120，N4=160，N5=200.

圖6為不同滑動(dòng)窗口的IAE濾波器估計(jì)的量測(cè)噪聲協(xié)方差.由圖可見(jiàn)，在測(cè)線上進(jìn)行飛行測(cè)量時(shí)，飛機(jī)的運(yùn)動(dòng)狀態(tài)平穩(wěn)，此時(shí)量測(cè)噪聲協(xié)方差估計(jì)曲線變化平緩；當(dāng)飛機(jī)由一條測(cè)線轉(zhuǎn)飛向另一條測(cè)線時(shí)(A、B、C、D、E指示的部分)，飛機(jī)的運(yùn)動(dòng)狀態(tài)變化劇烈，量測(cè)噪聲協(xié)方差估計(jì)曲線出現(xiàn)較大的跳動(dòng)，而且曲線的跳動(dòng)頻率和幅度隨著滑動(dòng)窗口寬度的增加而減小，即滑動(dòng)窗口越寬，IAE對(duì)量測(cè)噪聲的估計(jì)越平穩(wěn)，反之，IAE估計(jì)的量測(cè)噪聲曲線跳變?cè)絼×?

統(tǒng)計(jì)不同滑動(dòng)窗口的IAE與改進(jìn)IAE的解算誤差，如表2所示.不難看出，在N=80時(shí)，IAE的解算輸出與解算參考標(biāo)準(zhǔn)的誤差最小，此時(shí)的滑動(dòng)采樣窗口寬度兼顧了系統(tǒng)的運(yùn)動(dòng)狀態(tài)變化和新息協(xié)方差估計(jì)的準(zhǔn)確性；而在其他采樣窗口寬度條件下，IAE的解算結(jié)果誤差均高于N=80的情況.這也印證了前述的滑動(dòng)窗口寬度選取的不宜太大也不宜太小.

注意到利用改進(jìn)IAE解算的重力異常與參考標(biāo)準(zhǔn)的均方誤差在1.021047 mGal，解算精度劣于N=80時(shí)IAE的解算結(jié)果，但優(yōu)于其他的滑動(dòng)窗口寬度的情況.因此，改進(jìn)的IAE可以減小由于滑動(dòng)采樣窗口寬度選取不當(dāng)引起的解算誤差.

圖6　IAE估計(jì)R變化曲線(N=40、80、120、160、200)Fig.6　Estimated R by IAE (N=40,80,120,160,200)

濾波方法比較點(diǎn)數(shù)最大誤差(mGal)最小誤差(mGal)平均誤差(mGal)均方誤差(mGal)改進(jìn)IAE195184.660951-3.9856210.7810011.021047IAE(N=40)195186.116302-3.9087150.8110201.069920IAE(N=80)195183.872649-3.0882630.6840910.889861IAE(N=120)195184.070119-3.6601320.8671471.063740IAE(N=160)195184.310235-3.9103800.9236401.104894IAE(N=200)195184.645537-3.6042080.9348591.102127

5結(jié)語(yǔ)

本文利用新息自適應(yīng)卡爾曼濾波器及改進(jìn)的新息自適應(yīng)卡爾曼濾波器實(shí)現(xiàn)航空重力異常解算，并得出如下結(jié)論：

(1) 新息自適應(yīng)卡爾曼濾波器的應(yīng)用，可以降低量測(cè)噪聲協(xié)方差選取不準(zhǔn)確造成的解算誤差，減弱了解算精度對(duì)量測(cè)噪聲先驗(yàn)統(tǒng)計(jì)信息的依賴性；

(2) 對(duì)新息自適應(yīng)卡爾曼濾波器進(jìn)行改進(jìn)，即在新息自適應(yīng)卡爾曼濾波器中設(shè)計(jì)多個(gè)不同采樣窗口的新息協(xié)方差估計(jì)器，實(shí)現(xiàn)新息協(xié)方差估計(jì)器的滑動(dòng)采樣窗口優(yōu)化.解算結(jié)果表明，改進(jìn)的新息自適應(yīng)卡爾曼濾波的應(yīng)用能夠減小由滑動(dòng)采樣窗口選取不當(dāng)造成的解算誤差.

為獲得精度更高、更有普遍意義的航空重力異常解算方法，建議進(jìn)一步加強(qiáng)對(duì)自適應(yīng)卡爾曼濾波器設(shè)計(jì)方法的研究.同時(shí)，還應(yīng)繼續(xù)加強(qiáng)對(duì)航空重力測(cè)量系統(tǒng)的研究，進(jìn)一步完善航空重力異常解算的卡爾曼濾波模型，提高解算精度.

致謝感謝中國(guó)國(guó)土資源航空物探遙感中心的熊盛青、周堅(jiān)鑫、郭志宏和周錫華對(duì)本研究的大力支持和幫助.

References

Bian H W, Jin Z H, Wang J P, et al. 2006. The innovation-based estimation adaptive Kalman filter algorithm for INS/GPS integrated navigation system.JournalofShanghaiJiaotongUniversity(in Chinese), 40(6): 1000-1003, 1009.Bolotin Y V, Golovan A A, Parusniov N A. 2003. Mathematical models of airborne gravimetry systems with aided inertial platforms.MoscowUniversityMathematicsBulletin, 58(4):13-23.Bolotin Y V, Popelensky M Y. 2007. Accuracy analysis of airborne gravity when gravimeter parameters are identified in flight.JournalofMathematicalSciences, 146(3):5911-5519.

Bolotin Y V. 2009. Mathematics behind GTGRAV. Moscow: Moscow Lomonosov State University.

Bolotin Y V, Golovan A A. 2013. Methods of inertial gravimetry.MoscowUniversityMechanicsBulletin, 68(5): 117-125.

Cai S K, Wu M P, Zhang K D. 2010. A comparison of digital lowpass FIR-filters in airborne gravimetry.GeophysicalandGeochemicalExploration(in Chinese), 34(1):74-78.

Cai S K, Wu M P, Zhang K D, et al. 2013. The first airborne scalar gravimetry system based on SINS/DGPS in China.ScienceChinaEarthSciences, 56(12): 2198-2208.

Guo Z H, Luo F, An Z F. 2007. Experimental researches on FIR lowpass digital filters based on window functions of airborne gravity data.GeophysicalandGeochemicalExploration(in Chinese), 31(6):568-571, 576.

Guo Z H, Xiong S Q, Zhou J X, et al. 2008. The research on quality evaluation method of test repeat lines in airborne gravity survey.ChineseJournalofGeophysics(in Chinese), 51(5):1538-1543.

Guo Z H, Luo F, Wang M, et al. 2011. The design and experiment of IIR lowpass digital filters for airborne gravity data.ChineseJournalofGeophysics(in Chinese), 54(8):2148-2153.

Li R. 2014. Research of gravity anomaly solution in airborne gravimetry based on Kalman filtering [Master′s thesis] (in Chinese). Beijing: China University of Geosciences Beijing.

Luo F, Guo Z H, Luo Y, et al. 2012. Experimental researches on FIR lowpass filter based on equiripple.GeophysicalandGeochemicalExploration(in Chinese), 36(5): 856-860.

Mohamed A H. 1999. Optimizing the estimation procedure in INS/GPS integration for kinematic applications [Ph. D. thesis]. Calgary: Canada Calgary University.

Mohamed A H, Schwarz K P. 1999. Adaptive Kalman Filtering for INS/GPS.JournalofGeodesy, 73(4): 193-203.

Pan Y B. 2012. The Kalman filter method of free air anomaly solution in airborne gravimetry [Master′s thesis] (in Chinese). Beijing: China University of Geosciences Beijing.

Qin Y Y, Zhang H Y, Wang S H. 2012. Kalman Filter and Principle of Integrated Navigation 2nded (in Chinese). Xi′an: Northwestern Polytechnical University Press.Sun Z M. 2004. Theory, methods and applications of airborne gravimetry [Ph. D. thesis] (in Chinese). Zhengzhou: PLA Information Engineering University.

Tan F J, Xu J N, Li A, et al. 2009. Innovation-based adaptive Kalman filter for accelerometer signal de-noising.JournalofDataAcquisitionandProcessing(in Chinese), 24(2): 227-231.Wang J B, Xiong S Q, Zhou X H, et al. 2009. The advances in the study of the airborne gravimetry system.GeophysicalandGeochemicalExploration(in Chinese), 33(4): 368-373.

Wang J B, Xiong S Q, Guo Z H, et al. 2010. Estimation of the vertical acceleration for the airborne gravimetry using Kalman smoothing.ProgressinGeophysics(in Chinese), 25(3): 968-974, doi: 10.3969/j.issn.1004-2903.2010.03.035.

Wang J B. 2010. Methodologies and technology of data processing for airborne gravimetry [Ph. D. thesis] (in Chinese). Beijing: China University of Geosciences Beijing.

Wang J B, Xiong S Q, Guo Z H, et al. 2011. Research on methods for estimating airborne gravity anomalies.GeophysicalandGeochemicalExploration(in Chinese), 35(4): 493-498.

Xiong S Q. 2009a. The present situation and development of airborne gravity and magnetic survey techniques in China.ProgressinGeophysics(in Chinese), 24(1): 113-117.

Xiong S Q. 2009b. The strategic consideration of the development of China′s airborne geophysical technology.GeologyinChina(in Chinese), 36(6):1366-1374.

Xiong S Q, Zhou X H, Guo Z H, et al. 2010. Theory, Method and Application of Airborne Gravity Prospecting (in Chinese).Beijing: Geology Publishing.

Yue X K, Yuan J P. 2005. An adaptive Kalman filtering algorithm based on maximum-likelihood criterion.JournalofNorthwesternPolytechnicalUniversity(in Chinese), 23(4): 469-474.Zhang C D. 2005. On the subject of airborne gravimetry and airborne gravity gradiometry.ChineseJournalofEngineeringGeophysics(in Chinese), 2(4):282-291.Zheng W, Zhang G B, Li R. 2014. Application of Kalman filtering and smoothing for airborne gravimetry.AppliedMechanicsandMaterials, 548-549:1192-1195.

Zhou W, Cai T J. 2013. Study on filtering methods of airborne gravity.AppliedMechanicsandMaterials, 333-335:516-521.

附中文參考文獻(xiàn)

卞鴻巍, 金志華, 王俊璞等. 2006. 組合導(dǎo)航系統(tǒng)新息自適應(yīng)卡爾曼濾波算法. 上海交通大學(xué)學(xué)報(bào), 40(6): 1000-1003, 1009.

蔡劭琨, 吳美平, 張開(kāi)東. 2010. 航空重力測(cè)量中FIR低通濾波器的比較. 物探與化探, 34(1): 74-78.

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

郭志宏, 熊盛青, 周堅(jiān)鑫等. 2008. 航空重力重復(fù)線測(cè)試數(shù)據(jù)質(zhì)量評(píng)價(jià)方法研究.地球物理學(xué)報(bào), 51(5):1538-1543.

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

李瑞. 2014. 基于卡爾曼濾波的航空重力異常解算的研究[碩士論文]. 北京: 中國(guó)地質(zhì)大學(xué)(北京).

羅峰, 郭志宏, 駱遙等. 2012. 航空重力數(shù)據(jù)的等波紋FIR低通濾波試驗(yàn). 物探與化探, 36(5): 856-860.

潘炎冰. 2012. 航重測(cè)量中解算自由空氣異常的卡爾曼濾波方法[碩士論文]. 北京: 中國(guó)地質(zhì)大學(xué)(北京).

覃方君, 許江寧, 李安等. 2009. 基于新息自適應(yīng)卡爾曼濾波的加速度計(jì)信號(hào)降噪. 數(shù)據(jù)采集與處理, 24(2): 227-231.

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

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

王靜波, 熊盛青, 周錫華等. 2009. 航空重力測(cè)量系統(tǒng)研究進(jìn)展. 物探與化探, 33(4): 368-373.

王靜波, 熊盛青, 郭志宏等. 2010. 利用Kalman平滑技術(shù)估算航空重力測(cè)量中的載體垂直加速度. 地球物理學(xué)進(jìn)展, 25(3): 968-974, doi: 10.3969/j.issn.1004-2903.2010.03.035.

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

王靜波, 熊盛青, 郭志宏等. 2011. 航空重力異常估計(jì)方法研究. 物探與化探, 35(4): 493-498.

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

熊盛青. 2009b. 發(fā)展中國(guó)航空物探技術(shù)有關(guān)問(wèn)題的思考. 中國(guó)地質(zhì), 36(6): 1366-1374.

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

岳曉奎, 袁建平. 2005. 一種基于極大似然準(zhǔn)則的自適應(yīng)卡爾曼濾波算法. 西北工業(yè)大學(xué)學(xué)報(bào), 23(4): 469-474.

張昌達(dá). 2005. 航空重力測(cè)量和航空重力梯度測(cè)量問(wèn)題. 工程地球物理學(xué)報(bào), 2(4): 282-291.

(本文編輯胡素芳)

Application research on adaptive Kalman filtering for airborne gravity anomaly determination

ZHENG Wei1, 2, ZHANG Gui-Bin1*

1DepartmentofGeophysicsandInformationTechnology,ChinaUniversityofGeosciences,Beijing100083,China2DepartmentofInformationEngineering,HenanInstituteofScienceandTechnology,XinxiangHenan453003,China

AbstractBased on the basic principle of airborne gravimetry, Kalman filtering model for airborne gravity anomaly determination is established and innovation-based adaptive estimation Kalman filter (IAE) is applied to estimate the airborne gravity anomaly when measurement covariance is unknown. At the same time, improved IAE filter installing several innovation covariance estimators is designed to reduce the difficulty of sampling window width selection. It is can be seen from the experimental result that IAE is able to reduce gravity anomaly determination errors caused by lack of statistical information of measurement noise. Moreover, the improved IAE filter can decrease the rate of determination errors introduced by the improper selection of sampling window width, and the error of gravity anomaly determined by improved IAE is nearly 1 mGal.

KeywordsAirborne gravity anomaly determination； Adaptive Kalman filtering； Innovation； Covariance of measurement noise

基金項(xiàng)目國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)(2011AA060501)資助.

作者簡(jiǎn)介鄭崴，男，1984年生，博士，主要從事航空重力測(cè)量技術(shù)的研究. E-mail:13667060@qq.com *通訊作者張貴賓，男，1958年生，教授. 主要從事地球物理反演、重磁電勘探等研究. E-mail:gbzhang@cugb.edu.cn

doi:10.6038/cjg20160410 中圖分類號(hào)P223,P631

收稿日期2015-09-14，2015-12-24收修定稿

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

Zheng W, Zhang G B. 2016. Application research on adaptive Kalman filtering for airborne gravity anomaly determination.ChineseJ.Geophys. (in Chinese),59(4):1275-1283,doi:10.6038/cjg20160410.