李新宇，周召發(fā)，張志利，常振軍，郝詩文

火箭軍工程大學 兵器發(fā)射理論和技術(shù)國家重點學科實驗室，陜西 西安 710025

以飛機為載體的航空重力測量技術(shù)已成為高效測定地球重力場信息的主要手段[1]。在動態(tài)重力測量中，測量結(jié)果會受到發(fā)動機振動、載體動態(tài)性、空氣擾動等因素的干擾，測量噪聲廣泛分布，嚴重影響了重力測量的精度[2]。開展航空重力測量數(shù)據(jù)處理方法研究，提高數(shù)據(jù)處理的精度，具有十分重要的意義。

為了抑制噪聲對實現(xiàn)重力信息的高精度提取的影響，眾多學者對航空重力數(shù)據(jù)后處理進行了深入研究。頻域有限沖激響應（FIR）低通濾波方法在實際中被廣泛應用，可以有效消除高頻噪聲[3]。國外提出了具有特色的濾波方法，如R.R.B.Von Frese 等[4]提出了波數(shù)相關濾波器（WCF），通過比較測量數(shù)據(jù)與協(xié)同數(shù)據(jù)頻譜的相關程度來實現(xiàn)噪聲分離；B.A.Alberts 等[5]提出了頻域加權(quán)的方法，對測量數(shù)據(jù)的不同頻段賦予不同的權(quán)值，以此來進行降噪處理；Y.V.Bolotin 等[6]在時域處理方法的基礎上，將地球重力場的非均勻一致性建模為多狀態(tài)隱性馬爾可夫模型，對結(jié)果的細節(jié)有所提升。國內(nèi)方面，羅鋒等[7]提出使用Kalman平滑算法開展航空重力測量數(shù)據(jù)處理方法研究，在國內(nèi)首次解算出穩(wěn)定平臺式航空重力異常數(shù)據(jù)，且數(shù)據(jù)處理精度達到國際先進水平。參考文獻[8]將經(jīng)驗模態(tài)分解應用到航空重力測量數(shù)據(jù)處理中，對已經(jīng)過FIR 低通濾波的航空重力測量數(shù)據(jù)進行深化處理，可以很好地消除由飛機動態(tài)運動引起的動態(tài)誤差，進一步提高測量精度。針對一直存在的動態(tài)效應剩余影響，黃謨濤等[9]基于AIC信息量準則和互相關分析方法提出了一種針對各類動態(tài)效應剩余誤差的通用補償模型，內(nèi)符合精度從原來的9.35mGal提升到1.01mGal。為了削弱各類誤差源的影響，歐陽永忠等[10]提出了一種相關分析法和測線網(wǎng)平差兩階段綜合誤差補償方法，實測數(shù)據(jù)處理結(jié)果驗證了該方法的有效性和可靠性。

通過測線網(wǎng)平差求解系統(tǒng)誤差改正量，對各測線重力測量值進行補償是處理航空重力測量數(shù)據(jù)的重要方法[11]，但現(xiàn)階段平差方法往往僅適用于規(guī)則的測線網(wǎng)[12]。航空重力測量中常常采用重復測線的測量方式來檢驗航空重力儀動態(tài)測量的重復一致性，此時測線網(wǎng)平差的方法便不再適用[13]。蔡劭琨等[14]提出了一種航空重力測量重復測線系統(tǒng)誤差自調(diào)整的方法，也僅僅可以解決兩條重復測線的數(shù)據(jù)處理，怎樣對多條重復測線進行系統(tǒng)平差還需要繼續(xù)深入研究。

因此，本文提出了基于經(jīng)驗模態(tài)分解（EMD）和半系統(tǒng)誤差調(diào)整的航空重力測量重復測線數(shù)據(jù)處理方法。首先，推導了航空重力測量的基本原理，并詳細分析了測量噪聲的來源及誤差特性。其次，提出了基于EMD的誤差分離方法，對航空重力測量數(shù)據(jù)的深化處理進一步提高動態(tài)重力測量的精度。最重要的是，基于半系統(tǒng)誤差調(diào)整模型提出了適用于多條重復測線的系統(tǒng)平差方法，填補了該領域的研究空白。最后，通過對實測數(shù)據(jù)的處理，驗證了本文所提數(shù)據(jù)處理方法的精度和有效性。

1 航空重力測量數(shù)學模型

航空重力測量的兩個基本問題是如何保持重力敏感器的穩(wěn)定指向及如何從慣性加速度中分離出重力加速度[15]。平臺式重力儀基于跟蹤當?shù)氐乩硐档姆€(wěn)定平臺，可以隔絕載體角運動的影響，可以有效解決第一個問題。第二個問題的解決依賴于差分GPS，由其可以獲取高精度的位置和速度信息，用以計算載體垂向加速度以及重力相關的改正項，進而分離得到重力信息。

1.1 數(shù)學原理

式（1）所述的比力方程是動態(tài)重力測量和慣性導航系統(tǒng)共同的理論基礎

式中，δg為重力異常；f為重力敏感器測量得到的比力；h為載體高度信息，h?為載體垂向加速度；γ為正常重力矢量的垂向分量；ωie為地球自轉(zhuǎn)角速度；VE，VN分別為載體相對于地球的東向和北向速度；φ為載體緯度信息；RN，RM分別為地球參考橢球的卯酉圈和子午圈的曲率半徑；后三項合記為厄特弗斯改正項δaE。

除厄特弗斯改正項外，航空重力測量還涉及正常重力改正、空間改正、偏心改正、水平加速度改正及垂直加速度改正等多個改正項。

1.2 噪聲誤差源分析

以飛機為載體的航空動態(tài)重力測量數(shù)據(jù)受到大量具有不同頻譜特性且分布廣泛的噪聲影響，本節(jié)依據(jù)噪聲特性進行分類簡述。

1.2.1 高頻噪聲

在實際航空重力測量中，重力敏感器的輸出包含由飛機發(fā)動機以及空氣湍流造成的隨機振動加速度，這部分噪聲的強度是重力異常的數(shù)萬倍，屬于高頻噪聲。因此，需要利用低通濾波的方法進一步去除重力異常信息中的高頻噪聲。

1.2.2 未知特性噪聲

測量條件中各種隨機因素產(chǎn)生的偶然誤差具有偶然性和隨機性，由載體動態(tài)性、空氣擾動等因素產(chǎn)生的噪聲也無法明確其誤差特性?；诮?jīng)驗模態(tài)分解的數(shù)據(jù)深化處理是消除這些誤差的有效方法。

1.2.3 低頻誤差

受測量條件中某些特定因素如長周期傳感器誤差的系統(tǒng)性影響而產(chǎn)生的低頻誤差被稱為系統(tǒng)誤差，其與重力信息在同一頻段，無法通過低通濾波予以消除。系統(tǒng)誤差對觀測成果具有累積作用，因此在測量過程中常使用平差方法予以消除。

2 基于EMD的誤差分離

構(gòu)建補償模型是抑制重力測量動態(tài)效應影響的一種方法，但模型構(gòu)建的科學性及模型參數(shù)的準確性嚴重制約上述方法的有效性?；跀?shù)據(jù)處理的誤差分離是另一種有效的抑制方法。本節(jié)提出EMD 對已經(jīng)低通濾波過后的航空重力測量數(shù)據(jù)進行深化處理，削弱未知特性噪聲對重力測量結(jié)果的影響，進一步提高動態(tài)重力測量的精度。

2.1 EMD基本原理

EMD 方法認為，任何復雜的時間序列x(t)都可以分解為一組從高頻到低頻的本征模態(tài)函數(shù)（IMF）和一個殘余信號r(t)，這些本征模態(tài)函數(shù)可以很好地反映信號在任何時間局部的頻率特征[16]。分解結(jié)果如式（4）所示

每個IMF分量都滿足兩個條件：一是整個信號上極值點數(shù)目與零點數(shù)目相差為1；二是極大值與極小值分別構(gòu)成的包絡線均值為0。EMD的分解步驟為：（1）確定原序列x(t)中的極大值和極小值點，并采用三次樣條分別對極大值和極小值點進行插值，構(gòu)造上下包絡線xmax(t)和xmin(t)，可得上下包絡線的均值為m(t)；（2）計算原時間序列與包絡線均值序列的差值為h(t) =x(t)-m(t)；（3）若h(t)滿足IMF的兩個基本條件，則將h(t)設定為一個IMF分量，并求出剩余分量r(t) =x(t)-h(t)，繼續(xù)執(zhí)行步驟（4）；若不滿足，將h(t)替代x(t)返回執(zhí)行步驟（1）～（3）；（4）將r(t)作為新的序列重復執(zhí)行步驟（1）～（4）；循環(huán)執(zhí)行直到提取的最后一個IMF分量或者余項小于預先設定的閾值或已成為單調(diào)函數(shù)即可結(jié)束算法。

2.2 誤差分離

基于EMD對航空重力測量數(shù)據(jù)進行深化處理的誤差分離，首先對低通濾波后去除高頻噪聲的重力數(shù)據(jù)進行EMD分解，而后根據(jù)重力數(shù)據(jù)與IMF分量的相關性進行篩選并重構(gòu)得到去噪后的信號。誤差分離的流程如圖1所示。

圖1 基于EMD的誤差分離流程Fig.1 Error separation flowchart based on EMD

若IMF 分量與重力數(shù)據(jù)的相關系數(shù)r（t）> 0.7，則認定為強相關；若0.2

3 重復測線系統(tǒng)誤差調(diào)整

航空重力測量中各條重復測線在同一測點上的重力異常值理論上應該是相同的，但由于受到系統(tǒng)誤差的影響，實際測量值并不相等，存在測點不符值。本節(jié)基于測區(qū)半系統(tǒng)差調(diào)整的算法提出適用于重復測線的平差方法，對測點數(shù)據(jù)進行補償，抑制系統(tǒng)誤差的影響。

n條重復測線上各選取m個測點進行平差計算，重復測線上重力測量值分別為g1j，g2j，…，gnj，其中j表示某一測點，且j∈[1,m]。將上述測量值用以下模型表示。

式中，g0是j測點上的重力真值；εij和Δij分別表示測線1～n上的系統(tǒng)誤差和隨機誤差。對n條測線上的重力測量值計算均方誤差，其中測線i的均方誤差計算式為

選取均方誤差最小的測線上的測量值為重復測線的重力基準值gˉij。因此可得n條重復測線j測點上的重力不符值為

除去基準測線外的（n-1）條重復測線測點不符值的均值作為該測點的總體不符值，并以該值作為平差的指標，迭代計算至測點總體不符值代數(shù)和為0 時結(jié)束調(diào)整，并按照該方法遍歷重復測線上所有測點即可完成重復測線的系統(tǒng)平差處理。

4 實測數(shù)據(jù)分析

在某次航空重力測量中4條重復測線上獲得的重力異常結(jié)果如圖2所示，所示結(jié)果已通過相同的FIR低通濾波器將高頻噪聲濾除。由圖2 可以看出，重力異常數(shù)據(jù)低通濾波后相對平滑，高頻噪聲濾除效果較好，但仍然存在明顯的振蕩和測點不符值，計算得到內(nèi)符合精度為2.62mGal。

圖2 重復測線重力異常結(jié)果圖Fig.2 Gravity anomaly result map of repeated survey lines

4.1 EMD分解

前文已經(jīng)介紹了EMD分解的原理，下面首先對重復測線1 的重力異常結(jié)果進行EMD 分解，得到如圖3 所示的分解結(jié)果。進行相關性分析后選擇相關系數(shù)大于0.7 的IMF分量對重力異常信息進行重構(gòu)，得到EMD去噪前后重力異常結(jié)果對比圖，如圖4 所示，可知EMD 深化處理后重力異常結(jié)果中的振蕩得到了很好的分離。

圖3 reline1重力異常結(jié)果EMD分解結(jié)果圖Fig.3 EMD decomposition result diagram of reline1 gravity anomaly results

圖4 reline1重力異常數(shù)據(jù)深化處理前后對比Fig.4 Comparison between reline1 gravity anomaly data before and after deepening processing

使用同樣的方法對其他測線進行EMD分解并重構(gòu)，得到EMD深化處理后的重復測線重力異常測量結(jié)果，如圖5所示。重復測線數(shù)據(jù)EMD處理后消除了原有明顯的振蕩，計算可得內(nèi)符合精度為0.75mGal，但測點不符值依然存在，有待進一步平差處理。

圖5 EMD去噪后的重復測線重力異常結(jié)果圖Fig.5 Gravity anomaly results of repeated survey lines after EMD denoising

4.2 重復測線平差調(diào)整

在重復測線上等間隔選取756 個測點，按照前文設計的重復測線平差算法進行系統(tǒng)誤差調(diào)整。本文設計的系統(tǒng)平差算法精度與平差次數(shù)呈正相關，兩者的相關關系如圖6所示。由圖6 可知，迭代計算至第16 次后測點不符值精度已收斂，756 個測點不符值代數(shù)和以及標準差基本為0，此時認為平差結(jié)束。

圖6 平差算法精度與平差次數(shù)關系圖Fig.6 Relationship between adjustment algorithm accuracy and adjustment times

平差結(jié)束后的重復測線重力異常結(jié)果如圖7 所示，相較平差前的結(jié)果，單個測點的不符值仍然存在，似乎表示該算法的優(yōu)化作用不是特別明顯，但通過計算可知，此時重復測線的內(nèi)符合精度提升至0.28mGal，相較平差前的0.75mGal已有大幅提升。這是因為本文所提的重復測線平差方法是針對總體測線而言的，因此重復測線總體的內(nèi)符合精度提升較為明顯。

圖7 平差后重復測線重力異常結(jié)果圖Fig.7 Gravity anomaly results of repeated survey lines after adjustment

圖8給出了平差前后測點不符值的分布情況。顯而易見，平差結(jié)束后測點不符值誤差明顯減小，系統(tǒng)誤差得到有效補償。且測點不符值分布大致滿足標準正態(tài)分布，說明了各測點之間相互獨立。

圖8 平差前后測點不符值分布圖Fig.8 Distribution map of inconsistent values of measurement points before and after adjustment

4.3 結(jié)果分析

基于EMD分解有效分離了信號與誤差，本文提出的重復測線平差方法在總體上對系統(tǒng)誤差進行調(diào)整，進一步提高了總體測線的內(nèi)符合精度。重復測線內(nèi)符合精度統(tǒng)計見表1，由此可知本文的研究內(nèi)容顯著提升了數(shù)據(jù)處理的精度。

表1 重復測線內(nèi)符合精度統(tǒng)計Table 1 Statistical table of inner coincidence accuracy ofrepeated lines

5 結(jié)論

本文針對航空動態(tài)重力測量噪聲分布廣泛的問題，對重復測線的誤差分離展開研究，提出了基于EMD和半系統(tǒng)誤差調(diào)整的航空重力測量重復測線數(shù)據(jù)處理方法。實測數(shù)據(jù)的處理結(jié)果表明，EMD 的深化處理將內(nèi)符合精度從2.62mGal 提升至0.75mGal，有效分離出噪聲的影響；經(jīng)平差處理后抑制系統(tǒng)誤差的影響，內(nèi)符合精度進一步提升至0.28mGal，精度依次提升了71.37%和62.67%。本文研究內(nèi)容對進一步推動航空動態(tài)重力測量精度提升和工程應用化具有積極意義。