航空重力弱信號提取的卡爾曼濾波方法研究

王冠鑫，羅鋒，周錫華，閆方

(1.自然資源部 航空地球物理與遙感地質(zhì)重點實驗室，北京 100083; 2.中國自然資源航空物探遙感中心，北京 100083; 3. 北京自動化控制設備研究所,北京 100074)

0 引言

重力勘探是一種十分重要的地球物理勘探手段，主要包括航空重力測量、地面重力測量和海洋重力測量，當遇到復雜山地、沼澤、海洋等人工無法通過的情況或需要進行快速掃面測量時，通常采用航空重力測量對重力場數(shù)據(jù)進行采集。航空重力測量是以飛機為載體，機載重力儀、定位儀等設備，可快速、經(jīng)濟、綠色地獲得航空重力測量數(shù)據(jù)，在基礎地質(zhì)研究、油氣礦產(chǎn)資源勘查和國防建設等方面發(fā)揮了重要作用。但由于飛機自身的高頻振動、機體運動和氣流擾動等因素的影響，無論是捷聯(lián)式航空重力測量系統(tǒng)還是三軸穩(wěn)定平臺式(后文簡稱“平臺式”)航空重力測量系統(tǒng)獲得的重力儀垂向加速度計信號往往含有大量的擾動噪聲，而我們所要提取的重力異常信息十分微弱，信噪比高達幾千甚至上萬級，因此將航空重力異常稱為弱信號，借以體現(xiàn)其信號強度低和難以提取的數(shù)據(jù)特點。目前從原始測量數(shù)據(jù)中獲取重力異常信號(弱信號)成為一個關鍵技術難題，正嚴重制約著國內(nèi)外航空重力測量精度。因此，發(fā)展航空重力數(shù)據(jù)解算技術十分必要。

俄羅斯GT公司(gravimetric technology)生產(chǎn)的GT系列三軸穩(wěn)定平臺式航空重力測量系統(tǒng)和加拿大SGL公司(sander geophysics limited)的AIRGrav航空重力測量系統(tǒng)是當今世界上測量精度最高的系統(tǒng)，這些系統(tǒng)配有相應的重力數(shù)據(jù)處理軟件(“黑匣子”)，數(shù)據(jù)處理方法各具特色，不盡相同。其中，GT系列在航空重力數(shù)據(jù)處理中采用了卡爾曼濾波技術，解算精度高，尤其在抗擾動噪聲方面表現(xiàn)突出。

目前對于重力數(shù)據(jù)的信噪分離，國內(nèi)普遍采用FIR窗函數(shù)低通濾波[1-2]，主要根據(jù)處理人員的經(jīng)驗來判斷測量數(shù)據(jù)的質(zhì)量進而設定濾波器的截止頻率。雖然這種數(shù)據(jù)處理模式能夠從原始數(shù)據(jù)中提取出低頻段的有用信號，但由于航空重力場屬于位場，其信號分布于整個頻段，僅從頻率的角度對數(shù)據(jù)進行處理，勢必要舍棄高頻部分的有用信號，可能無法準確地解算航空重力異常。2010年起，國內(nèi)一些團隊[3-6]相繼開展了卡爾曼濾波技術在航空重力領域里的研究工作，雖然在理論研究和仿真測試中取得了一定的進展，但是目前仍沒有達到實際應用的需求。

通過對過往研究成果分析發(fā)現(xiàn)，常用于慣性系統(tǒng)解算的卡爾曼濾波算法是在時間域進行的，其在動態(tài)信號估計和預測方面，較傳統(tǒng)的頻域濾波“一刀切”的做法更具優(yōu)勢。而在航空重力弱信號的提取中，最重要的一部分工作就是將飛機運動過程中對重力加速度計造成的影響進行剝離，一旦針對某一測量系統(tǒng)建立了準確的濾波狀態(tài)方程且獲得了測量系統(tǒng)的參數(shù)，即可以對任意測區(qū)的數(shù)據(jù)進行解算，極大削弱了對數(shù)據(jù)處理人員經(jīng)驗的依賴程度，這充分說明了卡爾曼濾波在航空重力弱信號中的應用前景。

本文將圍繞卡爾曼濾波方法開展研究工作，在時域下對航空重力信號進行處理，無需確定信號的截止頻率，進而規(guī)避在數(shù)據(jù)精度和分辨率之間取舍的矛盾。卡爾曼濾波解算結(jié)果的質(zhì)量往往由兩方面決定：一是根據(jù)系統(tǒng)所列狀態(tài)方程的準確性；二是所選取的參數(shù)(先驗信息)與系統(tǒng)真實參數(shù)的偏離程度。而這兩個方面均需要對測量原理以及測量硬件特性有著充分的理解，這也是卡爾曼濾波技術始終沒有在我國航空重力測量領域完成實用化的主要原因。筆者將調(diào)整有確定控制的通用卡爾曼濾波公式，建立航空重力異常的數(shù)學模型，提出卡爾曼濾波狀態(tài)方程，解決重力信號與差分GNSS信號匹配、航空重力弱小信號提取的難題，旨在為后續(xù)針對各類新型國產(chǎn)航空重力測量系統(tǒng)的卡爾曼濾波狀態(tài)方程的構建提供理論依據(jù)和指導。研究者們只需在此基礎上，給出測量系統(tǒng)的先驗信息，即可將卡爾曼濾波運用到各類航空重力測量系統(tǒng)中。當然也可根據(jù)具體的測量系統(tǒng)對本文提出的狀態(tài)方程進行適應性調(diào)整，構建出更準確的卡爾曼濾波狀態(tài)方程。

1 有確定控制的離散卡爾曼濾波原理及修正

1.1 有確定控制的卡爾曼濾波公式及修正

設被估計狀態(tài)Xk在tk時刻受到系統(tǒng)噪聲序列Wk-1和確定性輸入項uk-1的驅(qū)動，則驅(qū)動機理的狀態(tài)方程和量測方程為

Xk=Φk,k-1Xk-1+Bk-1uk-1+Γk-1Wk-1，

(1)

Zk=HkXk+Vk，

(2)

式中：Φk,k-1、Bk-1為常數(shù)矩陣；uk-1為控制項；Γk-1為系統(tǒng)噪聲驅(qū)動陣；Hk為量測陣；Wk為系統(tǒng)激勵噪聲序列，其方差為Qk；Vk為量測噪聲序列，其方差為Rk；且Wk和Vk為互不相關白噪聲，其中假設Qk為非負定，Rk為正定[7]。

在航空重力實際測量中，其基本測量原理在該領域基本達成共識[8-9]。根據(jù)牛頓第二定律對測量載體的受力及飛行狀態(tài)進行分析，可構建航空重力異常的數(shù)學模型為

(3)

Xk=Φk,k-1Xk-1+Bkuk+ΓkWk。

(4)

實際上，重力異常近似于一個平穩(wěn)均勻的隨機過程，可以認為根據(jù)式(3)所構建的狀態(tài)方程的系統(tǒng)誤差為一個恒定的狀態(tài)[10]，則有

Xk=Φk,k-1Xk-1+Bkuk+ΓW，

(5)

式中：k=1,2,…,n。結(jié)合式(2)，即可構建航空重力離散卡爾曼濾波算法。

狀態(tài)一步預測：

(6)

狀態(tài)估計：

(7)

濾波增益：

(8)

一步預測均方誤差：

(9)

估計均方誤差：

(10)

1.2 最優(yōu)線性平滑

(11)

(12)

(13)

在此說明，本文所涉及到的所有卡爾曼濾波結(jié)果均為前向卡爾曼濾波和RTS平滑綜合處理的結(jié)果。

2 航空重力異常解算卡爾曼濾波狀態(tài)方程的提出

目前國內(nèi)所提出的卡爾曼濾波狀態(tài)方程均是針對特定系統(tǒng)所構建的，但是這類狀態(tài)方程[6]往往包含了材料疲勞程度、平臺安裝誤差角等某系統(tǒng)所特有的參數(shù)，針對性太強，移植到其他系統(tǒng)時會引起不適應性。為了提出一個適用于大多數(shù)航空重力測量系統(tǒng)的卡爾曼濾波狀態(tài)方程，為后續(xù)研究者提供參考，本文將忽略各系統(tǒng)的特性，僅根據(jù)測量載體在飛行過程中的運動狀態(tài)來構建狀態(tài)方程。

2.1 量測值的選取

由于航空重力測量系統(tǒng)存在差分GNSS解算的載體高度、速度及重力場等多個測量數(shù)據(jù)，在確定狀態(tài)方程前，首先要確定系統(tǒng)的量測值。圖1、2分別是某一架次實測改正后重力數(shù)據(jù)和載體高度數(shù)據(jù)的功率譜分析圖，從圖中可以看出重力數(shù)據(jù)的能量沒有明確的分布頻段，而高度數(shù)據(jù)的能量主要集中于低頻，后者更有利于R陣的選取。同時飛機運動變化產(chǎn)生的加速度信息會疊加在重力儀傳感器上，這就導致了重力數(shù)據(jù)的信噪比極低。由于差分GNSS解算的載體運動信息數(shù)據(jù)的信噪比相對較高，因此認為量測信息應來自于差分GNSS系統(tǒng)，而非機載重力儀系統(tǒng)，即將差分GNSS解算數(shù)據(jù)作為量測值。

圖1 改正后重力數(shù)據(jù)的功率譜分析Fig.1 Power spectrum analysis diagram of corrected gravity data

圖2 差分GNSS解算的載體高度數(shù)據(jù)功率譜分析Fig.2 Power spectrum analysis diagram of >differential GNSS height data

將航空重力異常數(shù)學模型式(3)變形可得用于描述航空重力測量系統(tǒng)的運動狀態(tài)方程

(14)

或

(15)

因此，觀測量可以為高度h或垂向速度vu。根據(jù)實際測試經(jīng)驗，將高度作為量測值進行解算時得到的重力異常噪聲含量更少。本文對某一架次差分GNSS解算的載體高度數(shù)據(jù)進行一次差分來獲得垂向速度數(shù)據(jù)(紅色)，并將其與該架次差分GNSS解算的載體垂向速度數(shù)據(jù)(藍色)進行對比(圖3)，從圖中可以明顯看出藍色曲線比紅色曲線的噪聲振幅要小，因此推薦優(yōu)先使用差分GNSS高度數(shù)據(jù)作為卡爾曼濾波的量測值。

圖3 差分GNSS解算的載體垂向速度(紅線)與高度一階導數(shù)(藍線)對比Fig.3 Comparison of vertical velocity (red line) and first derivative of height (blue line) calculated by differential GNSS

2.2 狀態(tài)方程的構建

為了構建連續(xù)卡爾曼濾波狀態(tài)方程，需對式(15)進行二次拆分，得到

(16)

重力異常是一個隨機過程[12]，可以用成型濾波器

(17)

(18)

2.3 先驗信息的選取問題

3 解算方法實際數(shù)據(jù)驗證

為了驗證所提方法的有效性和普適性，分別利用國產(chǎn)捷聯(lián)式航空重力測量系統(tǒng)、平臺式航空重力測量系統(tǒng)[13]的實際測量數(shù)據(jù)進行了卡爾曼濾波解算，提取出了高質(zhì)量的航空重力異常。

3.1 捷聯(lián)式航空重力測量數(shù)據(jù)解算應用

2016年12月，在某海域進行了國內(nèi)新研制的捷聯(lián)式航空重力測量系統(tǒng)重復線飛行測試。利用該次飛行數(shù)據(jù)對本文提出的卡爾曼濾波解算方法進行了測試，并與常用的FIR窗函數(shù)濾波器(截止頻率為100 s，后文統(tǒng)一稱為FIR100s濾波)的濾波結(jié)果進行了對比，利用重復線內(nèi)符合精度[14-16]對數(shù)據(jù)的處理質(zhì)量進行了評價(表1)，F(xiàn)IR100s濾波的重復線內(nèi)符合精度為0.659 mGal，卡爾曼濾波解算的重復線內(nèi)符合精度為0.605 mGal。圖4、圖5分別為捷聯(lián)式航空重力測量系統(tǒng)FIR100s濾波和卡爾曼濾波的解算結(jié)果。從圖中可以看出，兩者處理解算得到的航空重力異常曲線形態(tài)基本一致?？梢哉f明采用本文提出的卡爾曼濾波方法解算的航空重力異常結(jié)果準確可靠，且從重復線內(nèi)符合精度的角度來看，卡爾曼濾波解算方法略優(yōu)。

表1 捷聯(lián)式重力儀不同濾波方法重復線內(nèi)符合精度計算結(jié)果

圖4 FIR100s濾波的航空重力異常重復線內(nèi)符合精度Fig.4 The internal accord accuracy of repeat lines of FIR100s filter for airborne gravity anomaly

圖5 卡爾曼濾波的航空重力異常內(nèi)符合精度Fig.5 The internal accord accuracy of repeat lines of Kalman filter for airborne gravity anomaly

3.2 平臺式航空重力測量解算應用

2017年11月，在某海域進行了國產(chǎn)平臺式航空重力測量系統(tǒng)的飛行測試，完成了東西向重復線平飛及起伏飛行測試。利用該次飛行數(shù)據(jù)對本文提出的卡爾曼濾波解算方法進行了測試，并與FIR100s濾波結(jié)果進行了對比(表2)。所選用的數(shù)據(jù)為同一條測線的6次重復測量，包含了4條平飛(測線2701、2702、2703、2704)和2條機動飛行(“V”型起伏飛行(測線2705)、水平“S”型飛行(測線2706))。圖6、7分別為FIR100s濾波、卡爾曼濾波解算平飛測線的航空重力異常，前者的重復線內(nèi)符合精度為0.719 mGal，后者則為0.470 mGal。從圖中可以看出兩種濾波方法解算的航空重力異常曲線趨勢一致，但后者擾動噪聲明顯小于前者。

表2 平臺式重力儀不同濾波方法重復線內(nèi)符合精度計算結(jié)果

圖6 FIR100s濾波的航空重力異常內(nèi)符合精度(平飛)Fig.6 The internal accord accuracy of repeat lines of FIR100s filter for airborne gravity anomaly(in steady-state flight condition)

圖7 卡爾曼濾波的航空重力異常內(nèi)符合精度(平飛)Fig.7 The internal accord accuracy of repeat lines of Kalman filter for airborne gravity anomaly(in steady-state flight condition)

圖8、9分別為FIR100s濾波、卡爾曼濾波解算平飛和機動飛行測線的航空重力異常，前者的重復線內(nèi)符合精度為1.014 mGal，后者則為0.539 mGal。

圖8 FIR100s濾波的航空重力異常內(nèi)符合精度(平飛和機動)Fig.8 The internal accord accuracy of repeat lines of FIR100s filter for airborne gravity anomaly(in steady-state and maneuvering flight condition)

圖9 卡爾曼濾波的航空重力異常內(nèi)符合精度(平飛和機動)Fig.9 The internal accord accuracy of repeat lines of Kalman filter for airborne gravity anomaly(in steady-state and maneuvering flight condition)

從圖中可以看出，在機動飛行測試中，F(xiàn)IR100s的處理結(jié)果在測線上出現(xiàn)較大的起伏波動，而卡爾曼濾波則很好地去除了飛機機動飛行所造成的擾動影響，可以認為本文提出的卡爾曼濾波解算方法更能適應較為復雜的飛行情況，解算精度更高，異常更加可靠。

4 結(jié)論

1) 建立了航空重力異常的數(shù)學模型，對有確定控制的卡爾曼濾波公式進行了適應性調(diào)整，合理地將卡爾曼濾波方法應用到航空重力測量數(shù)據(jù)的解算中；并根據(jù)航空重力測量系統(tǒng)的測量原理，有針對性地提出了卡爾曼濾波狀態(tài)方程，并優(yōu)選差分GNSS高度數(shù)據(jù)作為卡爾曼濾波的量測值。

2) 利用實測數(shù)據(jù)對提出的卡爾曼濾波方法解算航空重力異常進行了應用測試，并與常規(guī)的FIR低通濾波(截止頻率為100 s)的結(jié)果進行了對比。測試結(jié)果表明：卡爾曼濾波不僅適用于捷聯(lián)式及平臺式航空重力測量數(shù)據(jù)的解算，解算出的航空重力異常精度高，而且在較為復雜的飛行情況下，比起FIR低通濾波方法，卡爾曼濾波方法解算的航空重力異常擾動少，精度高，異常可靠。

3) 研究和測試表明，從航空重力測量原理角度來建立卡爾曼濾波的狀態(tài)方程是簡單有效的，可以將其方便地應用到各類航空重力測量系統(tǒng)的數(shù)據(jù)處理解算上。

致謝：北京自動化控制設備研究所及國防科技大學智能科學學院為本次研究提供了試驗數(shù)據(jù)及部分技術支持，在此表示誠摯的謝意。