航空重力測量中垂直加速度改正的固定區(qū)間平滑算法研究*

鄭 崴,魏 輝,陳學(xué)鋒,閆雷兵

(1.河南工學(xué)院 電子信息工程學(xué)院,河南 新鄉(xiāng) 453003;2.新鄉(xiāng)市信號與信息處理重點實驗室,河南 新鄉(xiāng) 453003)

0 引言

航空重力測量是在飛行狀態(tài)下,利用航空重力測量系統(tǒng)進行重力測量的新型動態(tài)測量技術(shù)。測量過程中,航空重力儀會受到飛機發(fā)動機的振動、飛機垂向運動、氣流顛簸及飛機高度變化等造成的擾動加速度的干擾[1,2]。這些擾動加速度中又以飛機垂向運動對航空重力儀的干擾最大,其量級可達到104毫伽以上,遠大于百十毫伽重力異常信號[3,4]。并且,該擾動加速度具有一定的隨機性,無法采用嚴密的解析公式進行補償修正,因此,在航空重力測量數(shù)據(jù)處理階段,須采用濾波、平滑等多種信號處理技術(shù)來補償這部分擾動加速度干擾,這個過程即為航空重力數(shù)據(jù)處理的垂直加速度改正。目前,基于傅里葉變換的頻域數(shù)字濾波和基于系統(tǒng)模型的卡爾曼濾波是實現(xiàn)垂直加速度改正處理的主要方法[5-7]。我國研制的三軸穩(wěn)定平臺式航空重力測量系統(tǒng)的數(shù)據(jù)處理軟件即采用了卡爾曼濾波進行垂直加速度改正處理[8]?？柭鼮V波是在時域內(nèi)依據(jù)最優(yōu)估計理論實現(xiàn)全頻帶干擾的修正,依據(jù)狀態(tài)估計時刻所用到的量測信息情況,最優(yōu)估計可以分為預(yù)測、濾波和平滑三類。其中,濾波是利用當(dāng)前時刻以及此前時刻的所有量測信息對當(dāng)前狀態(tài)進行估計的算法。而平滑除了利用濾波所用的量測信息外還利用了當(dāng)前時刻以后的部分或所有時刻的量測信息[9]。因此,平滑是一種離線處理算法,能夠獲得優(yōu)于濾波的估計精度。常用的平滑算法主要有固定點平滑、固定滯后平滑和固定區(qū)間平滑。

結(jié)合航空重力測量垂直加速度改正事后處理的特點,平滑算法的應(yīng)用可以提高垂直加速度改正的精度。國內(nèi)外眾多學(xué)者開展了卡爾曼濾波在航空重力測量中的應(yīng)用研究。蔡體菁等人構(gòu)造了擴展卡爾曼濾波方程進行重力異常信號提取的方法[10];俄羅斯國立大學(xué)的Bolotin等人在慣性導(dǎo)航誤差補償模型的基礎(chǔ)上將待估計重力異常信號作為狀態(tài)向量構(gòu)建航空重力異常估計的模型,并通過實測數(shù)據(jù)驗證了構(gòu)建模型的有效性[11];張貴賓等構(gòu)建了基于重力異常模型的信息提取方法,并設(shè)計固定區(qū)間平滑器提高了重力異常估計精度[12,13]。盡管上述研究取得了一定的進展,但這些研究仍主要聚焦于濾波算法以及模型的構(gòu)建上,對于平滑算法研究較少。在具體的平滑算法、設(shè)計及處理方案等方面,對航空重力測量垂直加速度改正的平滑處理來說,還需進一步研究分析不同平滑算法在數(shù)據(jù)處理中的特性,從理論和應(yīng)用兩方面進行綜合考慮,選取合適的平滑算法。

本文針對航空重力測量數(shù)據(jù)處理中的垂直加速度改正處理,在經(jīng)典卡爾曼濾波算法的基礎(chǔ)上,結(jié)合工程實際條件,對兩種固定區(qū)間平滑算法——TFS(Tow-Filter-Smoothing)算法和RTS(Rung-Tung-Striebel)算法在重力異常提取中的應(yīng)用進行仿真試驗,分析兩種算法在航空重力測量數(shù)據(jù)處理中的特點和適用方案。

1 固定區(qū)間平滑算法原理

固定區(qū)間平滑是利用時間間隔內(nèi)所有量測值來估計系統(tǒng)在這個時間內(nèi)整個過程的狀態(tài)的算法[14]。TFS算法和RTS算法是兩種常用的固定區(qū)間平滑算法,其中TFS算法是先進行時間順序的卡爾曼濾波(前向濾波),然后再按照時間的逆序自后向前再次進行濾波(后向濾波),最后對兩次濾波結(jié)果進行數(shù)據(jù)融合實現(xiàn)平滑;RTS算法同樣也是先進行時間順序的卡爾曼濾波,但其后向濾波是在前向濾波的基礎(chǔ)上進行修正完成數(shù)據(jù)融合的。兩種平滑算法的流程圖如圖1所示。

(a)TFS算法流程

(b)RTS算法流程圖1 固定區(qū)間平滑算法流程圖

1.1 TFS算法

不難發(fā)現(xiàn),TFS算法實際上包含了兩個獨立的濾波過程——時間順序的前向濾波和時間逆序的后向濾波,之后再將兩個濾波結(jié)果進行數(shù)據(jù)融合,而數(shù)據(jù)融合則是以兩個濾波過程中估計的協(xié)方差為依據(jù)進行加權(quán)平均。其具體處理過程為:

第一部分,前向濾波。前向濾波本質(zhì)上就是標(biāo)準(zhǔn)卡爾曼濾波,對于一個線性離散系統(tǒng),其系統(tǒng)方程和量測方程可分別表示為

Xk=Φk,k-1Xk-1+Wk

(1)

Zk=HkXk+Vk

(2)

式中,Wk～N(0,Qk),Vk～N(0,Rk),則前向濾波按照k=1,2,…,N的順序依據(jù)卡爾曼濾波基本方程確定每個時刻的狀態(tài)估計[14]:

(3)

(4)

(5)

(6)

PFk=[I-KFkHk]PFk,k-1

(7)

第二部分,后向濾波。由于要保證前向和后向濾波的獨立性,后向濾波通常采用信息濾波算法按照k=N,N-1,…,1逆序進行狀態(tài)估計,有[14]

(8)

(9)

(10)

(11)

(12)

式中,IBk,IBk-1/k分別為后向估計協(xié)方差的逆和后向預(yù)測協(xié)方差的逆,即

(13)

Psk=(IFk+IBk)-1

(14)

(15)

1.2 RTS算法

(16)

(17)

(18)

可以看到,RTS算法是在后向處理過程中實現(xiàn)數(shù)據(jù)融合完成平滑的,相較于TFS平滑算法,其在后向濾波過程中同時實現(xiàn)數(shù)據(jù)融合。

2 數(shù)據(jù)仿真試驗

2.1 試驗設(shè)計

設(shè)待測區(qū)域有一半徑為400m、頂部埋深為100m的規(guī)則球體,且球體的剩余密度為0.15 g/cm3?，F(xiàn)有測量間隔為4m、長度為4000m的航空重力測量測線位于球體上空,球心位于測線2000m處,則可得重力異常理論值如圖2所示。

測量過程中,航空重力儀受到飛機隨機垂直運動產(chǎn)生的擾動加速度干擾,為獲得滿足精度要求的重力異常信息,需對測量結(jié)果進行濾波平滑處理。設(shè)垂直擾動加速度是均值為零、方差為20 mGal2的白噪聲,如圖3所示。不難看出,重力異常已經(jīng)完全被噪聲淹沒且噪聲的強度比所需信號高出一個數(shù)量級。

圖2 球體模型重力異常

圖3 重力異常與擾動噪聲干擾

2.2 平滑試驗與分析

現(xiàn)分別采用TFS平滑算法與RTS平滑算法對重力儀觀測結(jié)果進行濾波平滑處理。由于兩種平滑算法都基于卡爾曼濾波,因此需要建立系統(tǒng)方程和量測方程。設(shè)重力異常在采樣間隔的變化為ux,取其一次近似,有[13]

(19)

式中,x為測線位置,M為剩余質(zhì)量,D為球體模型埋深。采樣點處的重力異常為gx,采樣間隔為Δx,則可得系統(tǒng)方程

gx=gx-1+Δxux-1+qx-1

(20)

式中,qx-1是系統(tǒng)噪聲項,它是均值為零、方差為Qx-1的高斯白噪聲。同時,以重力儀的觀測方程作為量測方程,有

(21)

2.2.1 TFS平滑試驗

在構(gòu)建的系統(tǒng)方程和量測方程的基礎(chǔ)上,采用TFS算法進行平滑處理,前向濾波過程按照式(3)-(7)標(biāo)準(zhǔn)卡爾曼濾波方程進行處理并存儲所需的狀態(tài)估計和估計協(xié)方差,后向濾波按照式(8)-(13)進行處理并存儲后向濾波的狀態(tài)估計和估計協(xié)方差,之后再按照式(14)、(15)進行數(shù)據(jù)融合實現(xiàn)TFS平滑處理。平滑結(jié)果如圖4所示。

圖4 前向濾波、后向濾波和TFS平滑對比

由圖4可知,前向濾波結(jié)果由前向后逐漸收斂,而后向濾波結(jié)果由后向前逐漸收斂,但兩者相較于理論重力異常曲線都有較大的差別。而將前向濾波和后向濾波結(jié)果進行融合后所獲得的TFS平滑結(jié)果比較規(guī)則,也更接近于理論值。這說明,進行數(shù)據(jù)融合后所實現(xiàn)的平滑優(yōu)于兩次濾波的結(jié)果。

2.2.2 RTS平滑試驗

采用RTS算法進行平滑處理中,前向濾波同樣按時間順序依據(jù)卡爾曼濾波方程處理,并依次記錄所需的狀態(tài)預(yù)測、預(yù)測協(xié)方差、狀態(tài)估計及估計協(xié)方差。再按照式(16)—(18)實現(xiàn)平滑處理,平滑結(jié)果如圖5所示。不難看出,RTS算法同樣可以消除在濾波初始階段由于濾波器未收斂帶來的估計誤差,也使得經(jīng)過垂直加速度改正后的重力異常曲線更加平滑。

圖5 前向濾波和RTS平滑對比

2.2.3 TFS平滑與RTS平滑對比

將TFS算法平滑結(jié)果、RTS算法平滑結(jié)果分別同重力異常理論值進行對比,結(jié)果如圖6所示。不難看出,經(jīng)過TFS算法和RTS算法平滑處理后,都可以獲得與理論值形狀一致的重力異常結(jié)果,附加在測量結(jié)果的噪聲被壓制,平滑后的重力異常曲線形狀與理論值基本一致。

分別統(tǒng)計TFS算法與RTS算法平滑處理的誤差,結(jié)果如表1所示。TFS算法與RTS算法都可以較好地壓制擾動噪聲,兩者性能基本相當(dāng),TFS算法獲得重力異常均方誤差為0.006 mGal,RTS算法均方誤差為0.005 mGal,二者在精確度上略有差別。同時也應(yīng)注意到,TFS算法平滑處理的數(shù)據(jù)融合階段需要一直求解系統(tǒng)模型的逆矩陣,因此存在計算量大、逆向模型不易求解的問題。

圖6 TFS平滑與RTS平滑對比

表1 TFS平滑與RTS平滑處理重力異常誤差對比

3 結(jié)論

本文對固定區(qū)間平滑算法的基本原理進行了研究分析,并著重研究了TFS算法和RTS算法的原理,在此基礎(chǔ)上,將它們應(yīng)用于航空重力測量的模型試驗,并分析了兩種平滑算法的性能。

由試驗結(jié)果可知,在航空重力測量垂直加速度改正處理中,由于TFS算法和RTS算法都采用了更多的量測信息,可以獲得相較于濾波結(jié)果更精確的重力異常估計結(jié)果,提高了垂直加速度改正的精度。在對TFS算法和RTS算法的對比中,發(fā)現(xiàn)它們都可以對擾動噪聲進行較好的壓制,兩種算法的精度基本相當(dāng),RTS算法的精度略微優(yōu)于TFS算法。從性能、計算處理復(fù)雜程度等方面綜合考慮,在航空重力測量數(shù)據(jù)處理階段,可優(yōu)先采用RTS算法進行平滑處理。