蔡體菁，顏 穎，王新宇

(1.東南大學(xué) 儀器科學(xué)與工程學(xué)院，江蘇 南京 210096；2.中國電子科技集團(tuán)公司第二十六研究所，重慶400060)

0 引言

近十幾年來航空重力測量技術(shù)發(fā)展迅速，國際上不斷出現(xiàn)新一代航空重力儀，如GT-2A[1]。國內(nèi)也有多家科研單位研制出二軸、三軸穩(wěn)定平臺(tái)海空重力儀、捷聯(lián)式航空海洋重力儀[2-3]。在航空重力測量數(shù)據(jù)處理方面，?？罩亓xGT系列均采用了卡爾曼濾波技術(shù)，但具體算法未公開報(bào)道。2015年，東南大學(xué)給出了處理平臺(tái)式重力儀測量數(shù)據(jù)的平滑正、反卡爾曼濾波法[4]。近年來，一些作者對重力數(shù)據(jù)處理進(jìn)行了研究與討論[5-7]。本文針對國產(chǎn)捷聯(lián)式航空重力測量系統(tǒng)，根據(jù)馬爾可夫過程的重力異常模型，給出基于該模型的平滑綜合卡爾曼濾波法，它把速度作為觀測量，重力異常及其導(dǎo)數(shù)作為狀態(tài)變量，進(jìn)行組合導(dǎo)航和自由空間重力異常計(jì)算，并把計(jì)算結(jié)果與同機(jī)工作的GT-2A航空重力儀重力測量結(jié)果進(jìn)行比較，兩者測量誤差小于0.85 mGal。

1 平滑綜合的正反卡爾曼濾波

卡爾曼濾波是一種最小化方差估計(jì)的遞推數(shù)據(jù)處理算法，它通過前一時(shí)刻的狀態(tài)估計(jì)值和當(dāng)前時(shí)刻的觀測值，遞推出當(dāng)前時(shí)刻的狀態(tài)估計(jì)值。平滑是利用當(dāng)前時(shí)刻前、后的狀態(tài)值計(jì)算出當(dāng)前時(shí)刻的狀態(tài)值。常用的平滑方法有固定區(qū)間平滑、固定點(diǎn)平滑和固定滯后平滑3種。本文采用固定區(qū)間平滑方法。

隨機(jī)線性離散系統(tǒng)的狀態(tài)方程和觀測方程可寫為

(1)

式中：Xi為ti時(shí)刻的狀態(tài)向量；Φi+1，i為ti～ti+1時(shí)刻的狀態(tài)轉(zhuǎn)移矩陣；Bi為系統(tǒng)輸入矩陣；ui為系統(tǒng)輸入向量；Gi為系統(tǒng)噪聲輸入矩陣；wi為系統(tǒng)噪聲向量；Zi為ti時(shí)刻的觀測向量；Hi為觀測矩陣；vi為觀測噪聲向量。wi和vi均是零均值高斯白噪聲且相互獨(dú)立，其噪聲強(qiáng)度分別為Qi和Ri。

式(1)的正向卡爾曼濾波為

(2)

初始條件：

(3)

P0=P(t0)

(4)

式(1)的反向卡爾曼濾波為

(5)

初始條件：

(6)

(7)

式(1)的平滑綜合濾波為

(8)

2 GPS/捷聯(lián)慣性組合導(dǎo)航卡爾曼濾波方程

由于捷聯(lián)式航空重力儀和運(yùn)載體固連，其3個(gè)加速度計(jì)輸入軸方向會(huì)隨測量載體運(yùn)動(dòng)而不斷發(fā)生變化，因此，要獲得加速度計(jì)比力的垂向分量，就必須知道每一時(shí)刻運(yùn)載體的姿態(tài)，即載體坐標(biāo)系相對導(dǎo)航坐標(biāo)系的航向角、俯仰角和橫滾角。對于重力數(shù)據(jù)后處理而言，利用捷聯(lián)式航空重力儀和GPS基站的信息，應(yīng)用卡爾曼濾波進(jìn)行GPS/捷聯(lián)慣性組合導(dǎo)航計(jì)算是獲得高精度運(yùn)載體姿態(tài)最有效的方法。

GPS/捷聯(lián)慣性組合導(dǎo)航系統(tǒng)卡爾曼濾波狀態(tài)方程為

(9)

系統(tǒng)噪聲陣G=[ωgx,ωgy,ωgz,ωax,ωay,ωaz]T，其中ωgx、ωgy、ωgz分別為3個(gè)方向陀螺儀自身的白噪聲；ωax、ωay、ωaz分別為3個(gè)方向加速度計(jì)的一階馬氏過程驅(qū)動(dòng)白噪聲。式(9)的具體形式如下：

(11)

(11)

(12)

(13)

(14)

ε=εb+ωg

(15)

(16)

(17)

GPS/捷聯(lián)慣性組合導(dǎo)航系統(tǒng)卡爾曼濾波觀測方程為

z(t)=H(t)x(t)+v(t)

(18)

式中z=[LI-LG,λI-λG,hI-hG,VIE-VGE,VIN-VGN,VIU-VGU]T，下標(biāo)I為捷聯(lián)慣性導(dǎo)航計(jì)算值，G為GPS測量值。

3 重力異常卡爾曼濾波方程

捷聯(lián)式航空重力儀的輸出信息通常為載體坐標(biāo)系下的3個(gè)加速度、3個(gè)角速度、溫度和GPS移動(dòng)站的星歷、偽距、多普勒頻移、載波相位等原始數(shù)據(jù)。利用以上這些信息及GPS基站的原始信息，能夠計(jì)算得到測線上的自由空間重力異常。自由空間重力異常Δg可以通過對慣性導(dǎo)航比力方程的變形得到，其方程為

(19)

地球重力場重力異常模型可以用二階馬爾可夫過程來描述：

(20)

式(19)、(20)可寫為卡爾曼濾波狀態(tài)方程形式：

(21)

用差分GPS的速度信息作為卡爾曼濾波的觀測信息，觀測方程為

Vup=HX+ν

(14)

式中：Vup為垂向速度；H為測量矩陣；ν為測量噪聲。

4 計(jì)算結(jié)果

東南大學(xué)開展重力測量數(shù)據(jù)處理研究的捷聯(lián)式航空重力儀由高精度慣性測量單元、GPS接收機(jī)和數(shù)據(jù)采集器組成。慣性測量單元包含3個(gè)高精度單軸石英撓性擺式加速度計(jì)、3個(gè)高精度單軸激光陀螺儀、電源板、溫控板、計(jì)算機(jī)板、I/F板等電子線路，數(shù)據(jù)采集器采集慣性測量單元和GPS接收機(jī)輸出的信息，并保證信息精確同步。

根據(jù)捷聯(lián)式航空重力儀的測量數(shù)據(jù)，獲得測線上的重力異常，要進(jìn)行3步計(jì)算：

1)差分GPS位置、速度計(jì)算。

2)GPS/捷聯(lián)慣導(dǎo)組合導(dǎo)航解算。

3)重力異常解算。

航空重力測量數(shù)據(jù)處理流程如圖1所示。

圖1 航空重力測量數(shù)據(jù)處理流程圖

為了檢驗(yàn)平滑綜合的正、反卡爾曼濾波及重力異常模型的正確性和有效性，利用了捷聯(lián)式航空重力儀與俄羅斯的GT-2A航空重力儀同時(shí)飛行的測量數(shù)據(jù)。運(yùn)用上述計(jì)算方法得到自由空間重力異常，并把計(jì)算結(jié)果與GT-2A航空重力儀軟件處理的結(jié)果對比，對比結(jié)果如圖2所示。調(diào)整后，兩航空重力儀輸出的重力異常符合精度約為0.8 mGal。調(diào)整前的誤差主要是兩個(gè)航空重力儀的前?；鶞?zhǔn)點(diǎn)數(shù)據(jù)有差別，通常相差一個(gè)常數(shù)。由圖可見，在測線的細(xì)節(jié)上，GT-2A的計(jì)算結(jié)果更豐富。

圖2 兩航空重力儀輸出的重力異常結(jié)果比較

5 結(jié)束語

根據(jù)本文給出的GPS/捷聯(lián)慣性組合導(dǎo)航和把重力異常作為狀態(tài)量的卡爾曼濾波狀態(tài)方程，應(yīng)用平滑綜合的正、反卡爾曼濾波法處理捷聯(lián)式重力儀重力測量數(shù)據(jù)，能夠得到較滿意的重力異常值。