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

(1.河南工學(xué)院 電子信息工程學(xué)院,河南 新鄉(xiāng) 453003;2.新鄉(xiāng)市信號(hào)與信息處理重點(diǎn)實(shí)驗(yàn)室,河南 新鄉(xiāng) 453003)

0 引言

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

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

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

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

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

(a)TFS算法流程

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

1.1 TFS算法

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

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

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

(1)

Zk=HkXk+Vk

(2)

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

(3)

(4)

(5)

(6)

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

(7)

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

(8)

(9)

(10)

(11)

(12)

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

(13)

Psk=(IFk+IBk)-1

(14)

(15)

1.2 RTS算法

(16)

(17)

(18)

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

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

2.1 試驗(yàn)設(shè)計(jì)

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

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

圖2 球體模型重力異常

圖3 重力異常與擾動(dòng)噪聲干擾

2.2 平滑試驗(yàn)與分析

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

(19)

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

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

(20)

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

(21)

2.2.1 TFS平滑試驗(yàn)

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

圖4 前向?yàn)V波、后向?yàn)V波和TFS平滑對(duì)比

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

2.2.2 RTS平滑試驗(yàn)

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

圖5 前向?yàn)V波和RTS平滑對(duì)比

2.2.3 TFS平滑與RTS平滑對(duì)比

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

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

圖6 TFS平滑與RTS平滑對(duì)比

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

3 結(jié)論

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

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