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

鄭 崴,魏 輝,陳學鋒,閆雷兵

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

0 引言

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

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

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

1 固定區間平滑算法原理

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

(a)TFS算法流程

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

1.1 TFS算法

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

第一部分,前向濾波。前向濾波本質上就是標準卡爾曼濾波,對于一個線性離散系統,其系統方程和量測方程可分別表示為

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

(1)

Zk=HkXk+Vk

(2)

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

(3)

(4)

(5)

(6)

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

(7)

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

(8)

(9)

(10)

(11)

(12)

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

(13)

Psk=(IFk+IBk)-1

(14)

(15)

1.2 RTS算法

(16)

(17)

(18)

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

2 數據仿真試驗

2.1 試驗設計

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

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

圖2 球體模型重力異常

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

2.2 平滑試驗與分析

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

(19)

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

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

(20)

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

(21)

2.2.1 TFS平滑試驗

在構建的系統方程和量測方程的基礎上,采用TFS算法進行平滑處理,前向濾波過程按照式(3)-(7)標準卡爾曼濾波方程進行處理并存儲所需的狀態估計和估計協方差,后向濾波按照式(8)-(13)進行處理并存儲后向濾波的狀態估計和估計協方差,之后再按照式(14)、(15)進行數據融合實現TFS平滑處理。平滑結果如圖4所示。

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

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

2.2.2 RTS平滑試驗

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

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

2.2.3 TFS平滑與RTS平滑對比

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

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

圖6 TFS平滑與RTS平滑對比

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

3 結論

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

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