王冠鑫,羅鋒,周錫華,閆方
(1.自然資源部 航空地球物理與遙感地質(zhì)重點(diǎn)實(shí)驗(yàn)室,北京 100083; 2.中國自然資源航空物探遙感中心,北京 100083; 3. 北京自動(dòng)化控制設(shè)備研究所,北京 100074)
重力勘探是一種十分重要的地球物理勘探手段,主要包括航空重力測量、地面重力測量和海洋重力測量,當(dāng)遇到復(fù)雜山地、沼澤、海洋等人工無法通過的情況或需要進(jìn)行快速掃面測量時(shí),通常采用航空重力測量對重力場數(shù)據(jù)進(jìn)行采集。航空重力測量是以飛機(jī)為載體,機(jī)載重力儀、定位儀等設(shè)備,可快速、經(jīng)濟(jì)、綠色地獲得航空重力測量數(shù)據(jù),在基礎(chǔ)地質(zhì)研究、油氣礦產(chǎn)資源勘查和國防建設(shè)等方面發(fā)揮了重要作用。但由于飛機(jī)自身的高頻振動(dòng)、機(jī)體運(yùn)動(dòng)和氣流擾動(dòng)等因素的影響,無論是捷聯(lián)式航空重力測量系統(tǒng)還是三軸穩(wěn)定平臺式(后文簡稱“平臺式”)航空重力測量系統(tǒng)獲得的重力儀垂向加速度計(jì)信號往往含有大量的擾動(dòng)噪聲,而我們所要提取的重力異常信息十分微弱,信噪比高達(dá)幾千甚至上萬級,因此將航空重力異常稱為弱信號,借以體現(xiàn)其信號強(qiáng)度低和難以提取的數(shù)據(jù)特點(diǎn)。目前從原始測量數(shù)據(jù)中獲取重力異常信號(弱信號)成為一個(gè)關(guān)鍵技術(shù)難題,正嚴(yán)重制約著國內(nèi)外航空重力測量精度。因此,發(fā)展航空重力數(shù)據(jù)解算技術(shù)十分必要。
俄羅斯GT公司(gravimetric technology)生產(chǎn)的GT系列三軸穩(wěn)定平臺式航空重力測量系統(tǒng)和加拿大SGL公司(sander geophysics limited)的AIRGrav航空重力測量系統(tǒng)是當(dāng)今世界上測量精度最高的系統(tǒng),這些系統(tǒng)配有相應(yīng)的重力數(shù)據(jù)處理軟件(“黑匣子”),數(shù)據(jù)處理方法各具特色,不盡相同。其中,GT系列在航空重力數(shù)據(jù)處理中采用了卡爾曼濾波技術(shù),解算精度高,尤其在抗擾動(dòng)噪聲方面表現(xiàn)突出。
目前對于重力數(shù)據(jù)的信噪分離,國內(nèi)普遍采用FIR窗函數(shù)低通濾波[1-2],主要根據(jù)處理人員的經(jīng)驗(yàn)來判斷測量數(shù)據(jù)的質(zhì)量進(jìn)而設(shè)定濾波器的截止頻率。雖然這種數(shù)據(jù)處理模式能夠從原始數(shù)據(jù)中提取出低頻段的有用信號,但由于航空重力場屬于位場,其信號分布于整個(gè)頻段,僅從頻率的角度對數(shù)據(jù)進(jìn)行處理,勢必要舍棄高頻部分的有用信號,可能無法準(zhǔn)確地解算航空重力異常。2010年起,國內(nèi)一些團(tuán)隊(duì)[3-6]相繼開展了卡爾曼濾波技術(shù)在航空重力領(lǐng)域里的研究工作,雖然在理論研究和仿真測試中取得了一定的進(jìn)展,但是目前仍沒有達(dá)到實(shí)際應(yīng)用的需求。
通過對過往研究成果分析發(fā)現(xiàn),常用于慣性系統(tǒng)解算的卡爾曼濾波算法是在時(shí)間域進(jìn)行的,其在動(dòng)態(tài)信號估計(jì)和預(yù)測方面,較傳統(tǒng)的頻域?yàn)V波“一刀切”的做法更具優(yōu)勢。而在航空重力弱信號的提取中,最重要的一部分工作就是將飛機(jī)運(yùn)動(dòng)過程中對重力加速度計(jì)造成的影響進(jìn)行剝離,一旦針對某一測量系統(tǒng)建立了準(zhǔn)確的濾波狀態(tài)方程且獲得了測量系統(tǒng)的參數(shù),即可以對任意測區(qū)的數(shù)據(jù)進(jìn)行解算,極大削弱了對數(shù)據(jù)處理人員經(jīng)驗(yàn)的依賴程度,這充分說明了卡爾曼濾波在航空重力弱信號中的應(yīng)用前景。
本文將圍繞卡爾曼濾波方法開展研究工作,在時(shí)域下對航空重力信號進(jìn)行處理,無需確定信號的截止頻率,進(jìn)而規(guī)避在數(shù)據(jù)精度和分辨率之間取舍的矛盾。卡爾曼濾波解算結(jié)果的質(zhì)量往往由兩方面決定:一是根據(jù)系統(tǒng)所列狀態(tài)方程的準(zhǔn)確性;二是所選取的參數(shù)(先驗(yàn)信息)與系統(tǒng)真實(shí)參數(shù)的偏離程度。而這兩個(gè)方面均需要對測量原理以及測量硬件特性有著充分的理解,這也是卡爾曼濾波技術(shù)始終沒有在我國航空重力測量領(lǐng)域完成實(shí)用化的主要原因。筆者將調(diào)整有確定控制的通用卡爾曼濾波公式,建立航空重力異常的數(shù)學(xué)模型,提出卡爾曼濾波狀態(tài)方程,解決重力信號與差分GNSS信號匹配、航空重力弱小信號提取的難題,旨在為后續(xù)針對各類新型國產(chǎn)航空重力測量系統(tǒng)的卡爾曼濾波狀態(tài)方程的構(gòu)建提供理論依據(jù)和指導(dǎo)。研究者們只需在此基礎(chǔ)上,給出測量系統(tǒng)的先驗(yàn)信息,即可將卡爾曼濾波運(yùn)用到各類航空重力測量系統(tǒng)中。當(dāng)然也可根據(jù)具體的測量系統(tǒng)對本文提出的狀態(tài)方程進(jìn)行適應(yīng)性調(diào)整,構(gòu)建出更準(zhǔn)確的卡爾曼濾波狀態(tài)方程。
設(shè)被估計(jì)狀態(tài)Xk在tk時(shí)刻受到系統(tǒng)噪聲序列Wk-1和確定性輸入項(xiàng)uk-1的驅(qū)動(dòng),則驅(qū)動(dòng)機(jī)理的狀態(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為控制項(xiàng);Γk-1為系統(tǒng)噪聲驅(qū)動(dòng)陣;Hk為量測陣;Wk為系統(tǒng)激勵(lì)噪聲序列,其方差為Qk;Vk為量測噪聲序列,其方差為Rk;且Wk和Vk為互不相關(guān)白噪聲,其中假設(shè)Qk為非負(fù)定,Rk為正定[7]。
在航空重力實(shí)際測量中,其基本測量原理在該領(lǐng)域基本達(dá)成共識[8-9]。根據(jù)牛頓第二定律對測量載體的受力及飛行狀態(tài)進(jìn)行分析,可構(gòu)建航空重力異常的數(shù)學(xué)模型為
(3)
Xk=Φk,k-1Xk-1+Bkuk+ΓkWk。
(4)
實(shí)際上,重力異常近似于一個(gè)平穩(wěn)均勻的隨機(jī)過程,可以認(rèn)為根據(jù)式(3)所構(gòu)建的狀態(tài)方程的系統(tǒng)誤差為一個(gè)恒定的狀態(tài)[10],則有
Xk=Φk,k-1Xk-1+Bkuk+ΓW,
(5)
式中:k=1,2,…,n。結(jié)合式(2),即可構(gòu)建航空重力離散卡爾曼濾波算法。
狀態(tài)一步預(yù)測:
(6)
狀態(tài)估計(jì):
(7)
濾波增益:
(8)
一步預(yù)測均方誤差:
(9)
估計(jì)均方誤差:
(10)
(11)
(12)
(13)
在此說明,本文所涉及到的所有卡爾曼濾波結(jié)果均為前向卡爾曼濾波和RTS平滑綜合處理的結(jié)果。
目前國內(nèi)所提出的卡爾曼濾波狀態(tài)方程均是針對特定系統(tǒng)所構(gòu)建的,但是這類狀態(tài)方程[6]往往包含了材料疲勞程度、平臺安裝誤差角等某系統(tǒng)所特有的參數(shù),針對性太強(qiáng),移植到其他系統(tǒng)時(shí)會引起不適應(yīng)性。為了提出一個(gè)適用于大多數(shù)航空重力測量系統(tǒng)的卡爾曼濾波狀態(tài)方程,為后續(xù)研究者提供參考,本文將忽略各系統(tǒng)的特性,僅根據(jù)測量載體在飛行過程中的運(yùn)動(dòng)狀態(tài)來構(gòu)建狀態(tài)方程。
由于航空重力測量系統(tǒng)存在差分GNSS解算的載體高度、速度及重力場等多個(gè)測量數(shù)據(jù),在確定狀態(tài)方程前,首先要確定系統(tǒng)的量測值。圖1、2分別是某一架次實(shí)測改正后重力數(shù)據(jù)和載體高度數(shù)據(jù)的功率譜分析圖,從圖中可以看出重力數(shù)據(jù)的能量沒有明確的分布頻段,而高度數(shù)據(jù)的能量主要集中于低頻,后者更有利于R陣的選取。同時(shí)飛機(jī)運(yùn)動(dòng)變化產(chǎn)生的加速度信息會疊加在重力儀傳感器上,這就導(dǎo)致了重力數(shù)據(jù)的信噪比極低。由于差分GNSS解算的載體運(yùn)動(dòng)信息數(shù)據(jù)的信噪比相對較高,因此認(rèn)為量測信息應(yīng)來自于差分GNSS系統(tǒng),而非機(jī)載重力儀系統(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ù)學(xué)模型式(3)變形可得用于描述航空重力測量系統(tǒng)的運(yùn)動(dòng)狀態(tài)方程
(14)
或
(15)
因此,觀測量可以為高度h或垂向速度vu。根據(jù)實(shí)際測試經(jīng)驗(yàn),將高度作為量測值進(jìn)行解算時(shí)得到的重力異常噪聲含量更少。本文對某一架次差分GNSS解算的載體高度數(shù)據(jù)進(jìn)行一次差分來獲得垂向速度數(shù)據(jù)(紅色),并將其與該架次差分GNSS解算的載體垂向速度數(shù)據(jù)(藍(lán)色)進(jìn)行對比(圖3),從圖中可以明顯看出藍(lán)色曲線比紅色曲線的噪聲振幅要小,因此推薦優(yōu)先使用差分GNSS高度數(shù)據(jù)作為卡爾曼濾波的量測值。
圖3 差分GNSS解算的載體垂向速度(紅線)與高度一階導(dǎo)數(shù)(藍(lán)線)對比Fig.3 Comparison of vertical velocity (red line) and first derivative of height (blue line) calculated by differential GNSS
為了構(gòu)建連續(xù)卡爾曼濾波狀態(tài)方程,需對式(15)進(jìn)行二次拆分,得到
(16)
重力異常是一個(gè)隨機(jī)過程[12],可以用成型濾波器
(17)
(18)
為了驗(yàn)證所提方法的有效性和普適性,分別利用國產(chǎn)捷聯(lián)式航空重力測量系統(tǒng)、平臺式航空重力測量系統(tǒng)[13]的實(shí)際測量數(shù)據(jù)進(jìn)行了卡爾曼濾波解算,提取出了高質(zhì)量的航空重力異常。
2016年12月,在某海域進(jìn)行了國內(nèi)新研制的捷聯(lián)式航空重力測量系統(tǒng)重復(fù)線飛行測試。利用該次飛行數(shù)據(jù)對本文提出的卡爾曼濾波解算方法進(jìn)行了測試,并與常用的FIR窗函數(shù)濾波器(截止頻率為100 s,后文統(tǒng)一稱為FIR100s濾波)的濾波結(jié)果進(jìn)行了對比,利用重復(fù)線內(nèi)符合精度[14-16]對數(shù)據(jù)的處理質(zhì)量進(jìn)行了評價(jià)(表1),F(xiàn)IR100s濾波的重復(fù)線內(nèi)符合精度為0.659 mGal,卡爾曼濾波解算的重復(fù)線內(nèi)符合精度為0.605 mGal。圖4、圖5分別為捷聯(lián)式航空重力測量系統(tǒng)FIR100s濾波和卡爾曼濾波的解算結(jié)果。從圖中可以看出,兩者處理解算得到的航空重力異常曲線形態(tài)基本一致。可以說明采用本文提出的卡爾曼濾波方法解算的航空重力異常結(jié)果準(zhǔn)確可靠,且從重復(fù)線內(nèi)符合精度的角度來看,卡爾曼濾波解算方法略優(yōu)。
表1 捷聯(lián)式重力儀不同濾波方法重復(fù)線內(nèi)符合精度計(jì)算結(jié)果
圖4 FIR100s濾波的航空重力異常重復(fù)線內(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
2017年11月,在某海域進(jìn)行了國產(chǎn)平臺式航空重力測量系統(tǒng)的飛行測試,完成了東西向重復(fù)線平飛及起伏飛行測試。利用該次飛行數(shù)據(jù)對本文提出的卡爾曼濾波解算方法進(jìn)行了測試,并與FIR100s濾波結(jié)果進(jìn)行了對比(表2)。所選用的數(shù)據(jù)為同一條測線的6次重復(fù)測量,包含了4條平飛(測線2701、2702、2703、2704)和2條機(jī)動(dòng)飛行(“V”型起伏飛行(測線2705)、水平“S”型飛行(測線2706))。圖6、7分別為FIR100s濾波、卡爾曼濾波解算平飛測線的航空重力異常,前者的重復(fù)線內(nèi)符合精度為0.719 mGal,后者則為0.470 mGal。從圖中可以看出兩種濾波方法解算的航空重力異常曲線趨勢一致,但后者擾動(dòng)噪聲明顯小于前者。
表2 平臺式重力儀不同濾波方法重復(fù)線內(nèi)符合精度計(jì)算結(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濾波、卡爾曼濾波解算平飛和機(jī)動(dòng)飛行測線的航空重力異常,前者的重復(fù)線內(nèi)符合精度為1.014 mGal,后者則為0.539 mGal。
圖8 FIR100s濾波的航空重力異常內(nèi)符合精度(平飛和機(jī)動(dòng))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)符合精度(平飛和機(jī)動(dòng))Fig.9 The internal accord accuracy of repeat lines of Kalman filter for airborne gravity anomaly(in steady-state and maneuvering flight condition)
從圖中可以看出,在機(jī)動(dòng)飛行測試中,F(xiàn)IR100s的處理結(jié)果在測線上出現(xiàn)較大的起伏波動(dòng),而卡爾曼濾波則很好地去除了飛機(jī)機(jī)動(dòng)飛行所造成的擾動(dòng)影響,可以認(rèn)為本文提出的卡爾曼濾波解算方法更能適應(yīng)較為復(fù)雜的飛行情況,解算精度更高,異常更加可靠。
1) 建立了航空重力異常的數(shù)學(xué)模型,對有確定控制的卡爾曼濾波公式進(jìn)行了適應(yīng)性調(diào)整,合理地將卡爾曼濾波方法應(yīng)用到航空重力測量數(shù)據(jù)的解算中;并根據(jù)航空重力測量系統(tǒng)的測量原理,有針對性地提出了卡爾曼濾波狀態(tài)方程,并優(yōu)選差分GNSS高度數(shù)據(jù)作為卡爾曼濾波的量測值。
2) 利用實(shí)測數(shù)據(jù)對提出的卡爾曼濾波方法解算航空重力異常進(jìn)行了應(yīng)用測試,并與常規(guī)的FIR低通濾波(截止頻率為100 s)的結(jié)果進(jìn)行了對比。測試結(jié)果表明:卡爾曼濾波不僅適用于捷聯(lián)式及平臺式航空重力測量數(shù)據(jù)的解算,解算出的航空重力異常精度高,而且在較為復(fù)雜的飛行情況下,比起FIR低通濾波方法,卡爾曼濾波方法解算的航空重力異常擾動(dòng)少,精度高,異??煽俊?/p>
3) 研究和測試表明,從航空重力測量原理角度來建立卡爾曼濾波的狀態(tài)方程是簡單有效的,可以將其方便地應(yīng)用到各類航空重力測量系統(tǒng)的數(shù)據(jù)處理解算上。
致謝:北京自動(dòng)化控制設(shè)備研究所及國防科技大學(xué)智能科學(xué)學(xué)院為本次研究提供了試驗(yàn)數(shù)據(jù)及部分技術(shù)支持,在此表示誠摯的謝意。