修 睿，郭 剛，薛正兵，李東明，李海兵

（1.青島海洋科學(xué)與技術(shù)試點(diǎn)國(guó)家實(shí)驗(yàn)室,青島266237；2.北京航天控制儀器研究所,北京100039）

0 引言

航空重力測(cè)量是利用機(jī)載重力測(cè)量系統(tǒng)進(jìn)行的空中重力測(cè)量，是重力測(cè)量的重要手段。目前，國(guó)際上的航空重力測(cè)量系統(tǒng)主要為平臺(tái)式系統(tǒng)。近年來(lái)，國(guó)內(nèi)自主研制出了基于捷聯(lián)慣性技術(shù)的新型重力測(cè)量系統(tǒng)，其被稱之為捷聯(lián)式重力測(cè)量系統(tǒng)。捷聯(lián)式航空重力測(cè)量系統(tǒng)采用了捷聯(lián)慣性/GPS組合導(dǎo)航算法技術(shù)，可計(jì)算出重力傳感器敏感軸與當(dāng)?shù)氐乩硭矫娴牟淮怪倍炔⒂枰詳?shù)學(xué)校正，即可建立與物理平臺(tái)等價(jià)的數(shù)學(xué)平臺(tái)，具有體積小、質(zhì)量小、可靠性高的優(yōu)點(diǎn)。

在飛行作業(yè)過(guò)程中，由飛機(jī)機(jī)體高速運(yùn)動(dòng)帶來(lái)的高頻振動(dòng)不可避免地對(duì)重力傳感器產(chǎn)生了高頻干擾。為盡可能消除或減弱這些高頻噪聲的影響，需采用重力數(shù)據(jù)濾波方法對(duì)重力傳感器輸出進(jìn)行處理，從而獲得高精度的重力異常值。

目前，比較常用的濾波方法主要包括有限脈沖響應(yīng)FIR低通濾波、無(wú)限脈沖響應(yīng)IIR低通濾波、Kalman濾波、RC濾波、Gauss濾波等。本文采用FIR低通濾波、零相移濾波、標(biāo)準(zhǔn)Kalman濾波和正反Kalman濾波4種方法，結(jié)合實(shí)際飛行試驗(yàn)數(shù)據(jù)，對(duì)4種濾波方法的處理結(jié)果進(jìn)行了比對(duì)。同時(shí)，捷聯(lián)式重力測(cè)量系統(tǒng)的數(shù)學(xué)平臺(tái)精度也直接影響著重力傳感器敏感軸指向的校正精度，從而影響著最終的重力測(cè)量精度。本文對(duì)基于GPS輔助的SINS/GPS與基于DGPS輔助的SINS/DGPS 2種數(shù)學(xué)平臺(tái)計(jì)算方法獲得的重力測(cè)量數(shù)據(jù)結(jié)果進(jìn)行了比對(duì)。

1 系統(tǒng)組成及數(shù)據(jù)處理過(guò)程

1.1 捷聯(lián)式重力測(cè)量系統(tǒng)的測(cè)量原理

捷聯(lián)式?？罩亓x（SAG）將重力傳感器固聯(lián)集成于SINS/GPS組合導(dǎo)航系統(tǒng)，重力傳感器為高精度石英撓性加速度計(jì)。捷聯(lián)式?？罩亓x與顯控記錄裝置、UPS一起構(gòu)成了捷聯(lián)式?？罩亓y(cè)量系統(tǒng)。

重力傳感器敏感軸、捷聯(lián)慣性系統(tǒng)天向軸、載體天向軸三者平行，重力傳感器敏感軸的輸入信息包含了重力信息沿載體天向軸方向的分量及沿載體天向軸方向的運(yùn)動(dòng)加速度。通過(guò)SINS/GPS組合導(dǎo)航計(jì)算出載體姿態(tài)，即建立數(shù)學(xué)平臺(tái)，將重力傳感器輸出進(jìn)行指向校正，可獲得垂直于當(dāng)?shù)氐乩硭矫娴闹亓εc垂向運(yùn)動(dòng)加速度之和；再利用GPS系統(tǒng)測(cè)量信息估算出垂向運(yùn)動(dòng)加速度，即可獲得重力信息；再經(jīng)濾波處理濾除振動(dòng)等產(chǎn)生的干擾噪聲，以及厄缶改正、高度改正、正常重力改正等，最終可獲得自由空間重力異常值。

1.2 數(shù)據(jù)處理過(guò)程

捷聯(lián)式重力測(cè)量數(shù)據(jù)處理流程包括導(dǎo)航解算部分和重力異常提取部分，重力異常提取過(guò)程在導(dǎo)航解算的基礎(chǔ)上進(jìn)行。導(dǎo)航原始數(shù)據(jù)（陀螺、加速度計(jì)的原始脈沖）經(jīng)誤差補(bǔ)償和純慣性導(dǎo)航解算后，與GPS數(shù)據(jù)進(jìn)行組合導(dǎo)航和濾波處理，獲得純慣性導(dǎo)航的姿態(tài)誤差、速度誤差、位置誤差，修正得到東北天向的比力值和修正后的載體速度、位置信息。將載體位置信息進(jìn)行2次差分處理后，可得到載體的運(yùn)動(dòng)加速度。結(jié)合差分GPS速度、位置信息和3個(gè)比力信息，可以得到重力異常粗值。經(jīng)過(guò)濾波處理和各項(xiàng)改正后，可得到精確的自由空間重力異常值，其計(jì)算公式為[1]

式中，δg為自由空間重力異常值，gb為基準(zhǔn)點(diǎn)參考重力值，為比力和比力初值，為運(yùn)動(dòng)加速度改正，δaF為空間改正，δaE為厄缶改正，γ0為正常橢球面上的重力值。

自由空間重力異常提取由事后處理軟件完成，其流程如圖1所示。

圖1 自由空間重力異常提取處理流程圖Fig.1 Flowchart of free-space gravity anomaly extraction and processing

2 重力數(shù)據(jù)處理方法介紹

2.1 FIR低通濾波

目前，F(xiàn)IR濾波器主要涉及窗函數(shù)設(shè)計(jì)法、頻率抽樣設(shè)計(jì)法和Chebyshev逼近法，這里主要對(duì)常用的窗函數(shù)設(shè)計(jì)法進(jìn)行介紹。

有限脈沖響應(yīng)FIR數(shù)字濾波器要求用有限長(zhǎng)的單位沖擊響應(yīng)h（n）來(lái)逼近無(wú)限長(zhǎng)的理想濾波器的單位沖擊響應(yīng)hd（n），最常用和有效的方法就是用1個(gè)長(zhǎng)度為N的窗函數(shù)序列ω（n）來(lái)截取hd（n）的主要成分

實(shí)際上，窗函數(shù)設(shè)計(jì)FIR濾波器的核心是由給定的頻率特性通過(guò)加窗確定有限長(zhǎng)單位脈沖響應(yīng)序列h（n）。 工程中常見的窗函數(shù)包括以下幾種：矩形窗、三角窗、Hanning窗、Hamming窗、Blackman窗、Kaiser窗，其性能如表1所示[2-4]。

表1 常用窗函數(shù)的性能比較Table 1 Performance comparison of commonly used window functions

本文中用到的FIR濾波器采用Hanning窗，采樣頻率fs=1Hz，歸一化截止頻率fc=0.005Hz，階數(shù)N=600。

2.2 零相移濾波

零相移濾波是指信號(hào)序列在經(jīng)過(guò)濾波器濾波之后，信號(hào)序列相位不發(fā)生變化，系統(tǒng)函數(shù)的相位響應(yīng)為0，濾波后效果更平滑且信號(hào)無(wú)位移。零相移濾波算法流程如圖2所示，需進(jìn)行正向、反向2次濾波。

假設(shè)輸入信號(hào)序列為 [x（n）]，n=1，2，…，n0，[y（n）]為信號(hào)[x（n）]時(shí)域翻轉(zhuǎn)后的序列，將 2個(gè)信號(hào)分別進(jìn)行沿x軸、y軸的延拓后，進(jìn)行雙邊Z變換，得到

圖2 零相移算法原理圖Fig.2 Schematic diagram of zero phase shift algorithm

由算法流程，設(shè)傳遞函數(shù)為H（z），則濾波器輸出信號(hào)序列為

經(jīng)過(guò)二次濾波和翻轉(zhuǎn)輸出的信號(hào)序列為

本文中零相移濾波器的采樣頻率fs=1Hz，歸一化截止頻率fc=0.005Hz，階數(shù)N=600。

2.3 標(biāo)準(zhǔn)Kalman濾波

設(shè)tk時(shí)刻的被估計(jì)狀態(tài)Xk受系統(tǒng)噪聲序列Wk-1驅(qū)動(dòng)，其驅(qū)動(dòng)機(jī)理由下述狀態(tài)方程描述[5-7]

對(duì)Xk的量測(cè)滿足線性關(guān)系，量測(cè)方程為

式中，Φk，k-1為tk-1時(shí)刻至tk時(shí)刻的一步轉(zhuǎn)移陣，Γk-1為系統(tǒng)噪聲驅(qū)動(dòng)陣，Hk為量測(cè)陣，Vk為量測(cè)噪聲序列，Wk-1為系統(tǒng)激勵(lì)噪聲序列。

Kalman濾波方程為：

（1）狀態(tài)一步預(yù)測(cè)

（2）狀態(tài)估計(jì)

（3）濾波增益

（4）一步預(yù)測(cè)均方誤差

（5）估計(jì)均方誤差

式中，Pk為估計(jì)誤差方差陣，Kk為加權(quán)矩陣，Pk/k-1為一步預(yù)測(cè)均方誤差，Xk/k-1為狀態(tài)一步預(yù)測(cè)，Rk為量測(cè)噪聲方差陣。

2.4 正反Kalman濾波

正反Kalman濾波是連續(xù)進(jìn)行2次方向相反的Kalman濾波，并將正向、反向2次濾波結(jié)果進(jìn)行加權(quán)求和，能夠充分利用正向Kalman濾波的狀態(tài)估計(jì)量、一步狀態(tài)估計(jì)、誤差協(xié)方差矩陣等參數(shù)信息。該方法克服了標(biāo)準(zhǔn)離散Kalman濾波的精度缺陷問(wèn)題，優(yōu)化了濾波處理結(jié)果[8-10]。

正反Kalman濾波的正向過(guò)程即為標(biāo)準(zhǔn)離散型Kalman濾波。

反向過(guò)程的濾波方程為

預(yù)測(cè)過(guò)程方程為

加權(quán)過(guò)程方程為

正反Kalman狀態(tài)量是正向狀態(tài)估計(jì)量xk與反向狀態(tài)估計(jì)量xf（k/k-1）的加權(quán)求和，權(quán)重為Xzf（k）和Pzf（k）。

3 飛行試驗(yàn)結(jié)果比對(duì)

3.1 飛行試驗(yàn)概況

采用上述FIR低通濾波、零相移濾波、標(biāo)準(zhǔn)Kalman濾波、正反Kalman濾波4種濾波方法，對(duì)SAG某架次的飛行試驗(yàn)數(shù)據(jù)進(jìn)行了處理，并將處理結(jié)果進(jìn)行了精度比對(duì)。

在本次試驗(yàn)過(guò)程中，飛機(jī)平穩(wěn)飛行，飛行平均高度為800m，飛行平均速度為60m/s，飛行軌跡為8條東西重復(fù)測(cè)線、2條南北測(cè)線。測(cè)線結(jié)果如圖3所示，圖中的經(jīng)緯度進(jìn)行了人為偏移。

圖3 東西重復(fù)測(cè)線和南北測(cè)線結(jié)果Fig.3 Results of East-West repeat survey lines and North-South survey lines

3.2 基于SINS/GPS、SINS/DGPS組合導(dǎo)航的結(jié)果比對(duì)

分別采用GPS輔助導(dǎo)航和DGPS輔助導(dǎo)航進(jìn)行導(dǎo)航解算，2種組合導(dǎo)航解算得到的比力值、航姿、速度、位置在測(cè)線段的曲線分別如圖4～圖7所示，圖中的經(jīng)緯度進(jìn)行了人為偏移。其中，紅色曲線為GPS輔助導(dǎo)航的解算結(jié)果，藍(lán)色曲線為DGPS輔助導(dǎo)航的解算結(jié)果。對(duì)2種組合導(dǎo)航東西重復(fù)線上的比力、速度、航姿信息進(jìn)行均值和標(biāo)準(zhǔn)值的比對(duì)，其結(jié)果如表2所示。

圖4 2種組合導(dǎo)航解算測(cè)線段比力值的結(jié)果比對(duì)Fig.4 Comparison of specific force value between two kinds of integrated navigation solution survey lines

圖5 2種組合導(dǎo)航解算測(cè)線段航姿的結(jié)果比對(duì)Fig.5 Comparison of attitude angle between two kinds of integrated navigation solution survey lines

圖6 2種組合導(dǎo)航解算測(cè)線段速度的結(jié)果比對(duì)Fig.6 Comparison of velocity between two kinds of integrated navigation solution survey lines

圖7 2種組合導(dǎo)航解算測(cè)線段位置的結(jié)果比對(duì)Fig.7 Comparison of position between two kinds of integrated navigation solution survey lines

表2 2種組合導(dǎo)航測(cè)線段比力、姿態(tài)、速度的信息比對(duì)Table 2 Comparison of specific force，attitude angle and velocity between two kinds of integrated navigation solution survey lines

由圖4～圖7所示的2種組合導(dǎo)航的比力值、航姿、速度、位置結(jié)果的比對(duì)曲線可以看出，紅色曲線和藍(lán)色曲線均完全重合，即GPS輔助導(dǎo)航結(jié)果與DGPS輔助導(dǎo)航的結(jié)果基本完全一致。由表2的比對(duì)數(shù)據(jù)可以看出，各項(xiàng)信息的均值和標(biāo)準(zhǔn)值基本沒(méi)有差異，GPS輔助導(dǎo)航結(jié)果與DGPS輔助導(dǎo)航結(jié)果一致，采用單點(diǎn)GPS或差分GPS對(duì)于導(dǎo)航信息和比力值的坐標(biāo)分解沒(méi)有影響。

3.3 基于SINS/GPS、SINS/DGPS組合導(dǎo)航的重力數(shù)據(jù)處理結(jié)果比對(duì)

分別對(duì)基于SINS/GPS和SINS/DGPS組合導(dǎo)航的重力測(cè)量數(shù)據(jù)，采用FIR低通濾波、零相移濾波、標(biāo)準(zhǔn)Kalman濾波、正反Kalman濾波方法對(duì)測(cè)區(qū)重力異常粗值進(jìn)行濾波處理。其中，在重力數(shù)據(jù)處理階段的各項(xiàng)改正中代入的速度、位置、高度數(shù)據(jù)均為差分GPS數(shù)據(jù)。處理結(jié)果分別如圖8、圖9所示，藍(lán)色曲線為FIR低通濾波結(jié)果，紅色曲線為標(biāo)準(zhǔn)Kalman濾波結(jié)果，綠色曲線為正反Kalman濾波結(jié)果，黑色曲線為零相位濾波結(jié)果。

圖8 基于SINS/GPS的4種濾波方法處理重力異常結(jié)果比對(duì)Fig.8 Comparison of four filtering methods for gravity anomaly results based on SINS/GPS

圖9 基于SINS/DGPS的4種濾波方法處理重力異常結(jié)果比對(duì)Fig.9 Comparison of four filtering methods for gravity anomaly results based on SINS/DGPS

對(duì)8條東西重復(fù)測(cè)線上的重力異常值進(jìn)行濾波處理，結(jié)果比對(duì)如圖10、圖11所示。藍(lán)色曲線為FIR低通濾波結(jié)果，紅色曲線為標(biāo)準(zhǔn)Kalman濾波結(jié)果，綠色曲線為正反Kalman濾波結(jié)果，黑色曲線為零相移濾波結(jié)果。

為了更好地比較4種濾波方法及不同組合導(dǎo)航條件下的重力異常數(shù)據(jù)處理結(jié)果，表3列出了這8條東西重復(fù)測(cè)線在水平調(diào)整后的內(nèi)符合精度。其中，最后1列為所有測(cè)線的內(nèi)符合總精度。

圖10 基于SINS/GPS的8條東西重復(fù)測(cè)線的4種濾波方法比對(duì)Fig.10 Comparison of four filtering methods on eight kinds of East-West repeated survey lines based on SINS/GPS

圖11 基于SINS/DGPS的8條東西重復(fù)測(cè)線的4種濾波方法比對(duì)Fig.11 Comparison of four filtering methods on eight kinds of East-West repeated survey lines based on SINS/DGPS

表3 重復(fù)測(cè)線內(nèi)符合精度評(píng)估Table 3 Accuracy evaluation in repeated survey lines

在表3中，以每條重復(fù)測(cè)線的內(nèi)符合精度及所有重復(fù)測(cè)線的總內(nèi)符合精度進(jìn)行評(píng)判，分別對(duì)采用單點(diǎn)GPS進(jìn)行組合導(dǎo)航后的相關(guān)重力測(cè)量數(shù)據(jù)處理結(jié)果及采用差分GPS進(jìn)行組合導(dǎo)航后的相關(guān)重力測(cè)量數(shù)據(jù)處理結(jié)果進(jìn)行比對(duì)，并列出了采用4種濾波方法所得的重復(fù)測(cè)線內(nèi)符合精度，均為水平調(diào)整后的結(jié)果。

每條測(cè)線的均方根誤差及標(biāo)準(zhǔn)差計(jì)算方法如下

式中，δij為第j條重復(fù)公共段各點(diǎn)觀測(cè)值與該點(diǎn)各重復(fù)線觀測(cè)值的平均值之差，即

m為重復(fù)線數(shù)目，n為重復(fù)線公共段數(shù)據(jù)點(diǎn)數(shù)。

所有重復(fù)線測(cè)試數(shù)據(jù)的總內(nèi)符合精度計(jì)算公式為

從表3基于SINS/GPS組合導(dǎo)航的重力數(shù)據(jù)處理結(jié)果的最后一列，即所有重復(fù)測(cè)線的內(nèi)符合總精度可以看出：2種Kalman濾波方式相比，正反Kalman濾波的均方根誤差最小，為1mGal左右，標(biāo)準(zhǔn)Kalman濾波的均方根誤差最大，達(dá)到近1.38mGal，兩者相差0.38mGal；2種低通濾波相比較為接近，F(xiàn)IR低通濾波的均方根誤差為1.19mGal，零相移濾波的精度為1.08mGal，略優(yōu)于FIR低通濾波；而零相移濾波與正反Kalman濾波的結(jié)果相近，后者精度略高。

從表3基于SINS/DGPS組合導(dǎo)航的重力數(shù)據(jù)處理結(jié)果的最后1列數(shù)據(jù)可以看出：正反Kalman濾波效果最好，精度達(dá)到了0.98mGal；零相移濾波其次，達(dá)到了1.07mGal；FIR低通濾波的精度相對(duì)較差，為1.2mGal；標(biāo)準(zhǔn)Kalman濾波的精度結(jié)果為1.39mGal。

將兩者比對(duì)可以看出：采用差分GPS數(shù)據(jù)進(jìn)行組合導(dǎo)航處理之后得到的重力測(cè)量數(shù)據(jù)，相比于采用單點(diǎn)GPS數(shù)據(jù)進(jìn)行組合導(dǎo)航得到的重力測(cè)量數(shù)據(jù)，精度基本接近，相差約為0.01mGal。

從以上分析可以看出：正反Kalman濾波方法的處理精度最高，優(yōu)于2種低通濾波和單向標(biāo)準(zhǔn)Kalman濾波；正反Kalman濾波得到的重力異常曲線與另外3種曲線相比最為平滑，對(duì)高頻噪聲的去噪效果最好；零相移濾波與FIR低通濾波相比，其優(yōu)點(diǎn)在于可進(jìn)行正反2次低通濾波，對(duì)高頻噪聲的去噪效果更好，并抵消了正反相位移動(dòng)的偏差，相對(duì)而言結(jié)果更為準(zhǔn)確。同時(shí)，采用SINS/DGPS組合導(dǎo)航得到的重力測(cè)量數(shù)據(jù)，相比于采用SINS/GPS組合導(dǎo)航得到的重力測(cè)量數(shù)據(jù)，處理結(jié)果基本一致，這說(shuō)明采用單點(diǎn)GPS或差分GPS進(jìn)行組合導(dǎo)航，對(duì)導(dǎo)航解算結(jié)果和重力信息數(shù)據(jù)的處理結(jié)果并不影響。

4 結(jié)論

本文對(duì)FIR低通濾波、零相移濾波、標(biāo)準(zhǔn)Kalman濾波、正反Kalman濾波4種重力異常提取方法進(jìn)行了比對(duì)，并采用以上方法分別對(duì)基于SINS/GPS和SINS/DGPS組合導(dǎo)航的實(shí)際SAG飛行試驗(yàn)數(shù)據(jù)進(jìn)行了處理。重復(fù)測(cè)線的內(nèi)符合精度統(tǒng)計(jì)結(jié)果表明，正反Kalman濾波與零相移濾波的結(jié)果精度相對(duì)較高，分別達(dá)到了0.98mGal和1.07mGal，正反Kalman濾波略優(yōu)于零相移濾波；標(biāo)準(zhǔn)Kalman濾波的精度相對(duì)較差，精度與正反Kalman濾波相差近0.4mGal；FIR低通濾波的精度適中，與零相移濾波相比略有差距。本文驗(yàn)證了這4種濾波方法的實(shí)用性和優(yōu)劣性，以及重力異常模型和窗函數(shù)設(shè)計(jì)的合理性。同時(shí)，也驗(yàn)證了單點(diǎn)GPS輔助導(dǎo)航與DGPS輔助導(dǎo)航的結(jié)果相同，即采用何種組合導(dǎo)航方式不影響后續(xù)重力數(shù)據(jù)的處理精度。本次試驗(yàn)也從側(cè)面反映了SAG重力儀的測(cè)量精度，其內(nèi)符合精度可達(dá)0.98mGal。